T  o  L  l  I 1  v  qv; 


NUSC  Technical  Report  6659 


27  May  1982 


IFD:  An  Implicit  Finite-Difference 
Computer  Model  for  Solving 
the  Parabolic  Equation 


Ding  Lee 
George  Boteeas 

Surface  Ship  Sonar  Department 


Q_ 

O 

CO 

Li. 

Lil 


Naval  Underwater  Systems  Center 

Newport,  Rhode  Island/New  London,  Connecticut 


Approved  (or  public  rdoaco;  distribution  unlimited 


88  07  26  1^6 


Preface 


t 


This  report  was  prepared  under  NUSC  Project  No.  A65020,  Finite-Difference 
Solutions  to  Acoustic  Wave  Propagation,  Principal  Investigator  Dr.  D.  Lee,  Code 
3342.  This  work  was  supported  by  NAVMAT,  Program  Manager  CAPT  D.  F. 
Parrish,  Code  08L,  Program  Element  61152N,  Navy  Sub/Task  ZROOOOOIOI*^- 
House  Laboratory  Independent  Research. 

The  Technical  Reviewer  for  this  report  was  Dr.  A.  H.  Quazi,  Code  3331 . 


Reviewed  .ad  Approved:  27  May  1982 


ytiJ A 

/  W.  A.  Von  Winkle 
Associate  Technical  Director  for  Technology 


The  authors  of  this  report  are  located  at  the  New  London 
Laboratory,  Naval  Underwater  Systems  Center, 

New  London,  Connecticut  06320. 


REPORT  DOCUMENTATION  PAGE 


*  Ti  T'.E  'tnd  JuAdiUi 


/fp:  AN  IMPLICIT  FINITE-DIFFERENCE  COMPUTER  MODEL  FOR 
SOLVING  THE  PARABOLIC  EQUATION 


7  AUTNOA, 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


3-  AECIAlENT‘S  CATALOG  NUMAEA 


s  typs  or  aeaqat  »  acaioo  coveaeo 


«■  ASAFOAMINO  O AO.  acaoat  NUMEEA 


I.  CONTAACT  QA  OAANT  NUMEEAf*) 


Ding  Lee 
George  Botseas 


S.  AKATO AMINO  OAOANIZATION  NAME  ANO  AOOAESS 

Naval  Underwater  Systems  Center 

New  London  Laboratory,  New  London,  CT.  06320 


II  CONTAOLLINO  OFFICE  NAME  ANO  AOOAESS 

Naval  Material  Command,  Code  08L 
Washington,  D.C.  20362 

I  ion 


k.  MONITORING  AGENCY  name  6  AOOPES S(lt  different  from  Controlling  Office)  19.  SECURITY  CLASS,  (ol  thta  report) 

UNCLASSIFIED 


II.  NUMBER  OP  RAGES 


is*,  oeclassipication/oowngraoing 
schedule 


16.  DISTRIBUTION  STATEMENT  (Of  (hi a  Report) 


Approved  for  public  release;  distribution  unlimited. 


IT.  DISTRIBUTION  STATEMENT  (of  the  at  a  tract  entered  m  Block  20.  If  different  from  Report) 


»6.  SUPPLEMENT  ARY  NOTES 


19.  KEY  WOROS  (Continue  on  roworao  a  I  do  If  nocoooary  and  Identify  br 


Parabolic  wave  equation 
Imolicit  finite  difference 
Numerical  methods 
Propagation  loss 


I'..  ^r  acoustics 
Com,, model 
Range  independent 
Range  dependent 


20.  ABSTRACT  (Continue  on  reeereo  at  do  if  nocoeeary  and  identity  by  block  number) 

-'b  general  purpose  computer  model,  based  on  the  implementation  of  an 
implicit  finite-difference  scheme,  is  developed  for  the  solution  of  the 
parabolic  wave  equation.  This  model  can  be  used  to  predict  acoustic  propaga¬ 
tion  loss  in  both  range- dependent  and  range- independent  environments.  An 
important  feature  of  the  model  is  that  it  can  handle  arbitrary  surface  boundary 
conditions  and  an  irregular  bottom  with  arbitrary  bottom  boundary  conditions. 

In  the  event  that  the  bottom  boundary  conditions  cannot  be  expressed  >. 


1473  EOITION  OF  I  NOV  SS  IS  OMOLETl 

S/A  3  10  3-0  14*  «60  l 


on  wom* 

UW  I  J  AN  71 


20.  (Continued) 


mathematically,  the  model  has  the  capability  of  introducina  an  artifical 
absorbing  bottom  such  that,  with  the  aoprooriate  bottom  attenuation,  the 
problem  becomes  solvable.  Another  important  feature  of  the  model  is  that 
it  can  handle  horizontal  interfaces  of  layered  media.  The  model  is  easy 
to  use  and  easy  to  modify.  Numerical  test  examples  are  included  to  demon¬ 
strate  the  capabilities  of  the  model.  . 


TR  6659 


TABLE  OF  CONTENTS 


Page 


LIST  OF  ILLUSTRATIONS . ii 

1.  INTRODUCTION . 1-1 

2.  SOLUTION  BACKGROUND . 2-1 

2.1.  Theoretical  Derivation  .  2-1 

2.2.  Implicit  Finite-Difference  Formula  .  2-2 

2.3.  Boundary  Treatment  .  2-5 

2.4.  Interface  Treatment  .  2-6 

2.5.  Attenuation . 2-11 

3.  COMPUTER  IMPLEMENTATION  .  3-1 

3.1.  Desirable  Features  .  3-1 

3.2.  Structure  of  Program . 3-2 

3.3.  User's  Guide . 3-10 

4.  TEST  PROBLEMS . 4-1 

4.1.  Exact  Solution . 4-1 

4.2.  Range- Independent  Problems  .  4-8 

4.3.  Range- Dependent  Problems  .  4-11 

5.  CONCLUSIONS . 5-1 

6.  REFERENCES . 6-1 

APPENDIX  A — IFD  COMPUTER  LISTING  .  A-l 

APPENDIX  B— PLOT  PROGRAM  COMPUTER  LISTING  .  B-l 


Accession  For 

Nil?  GPAtl 
OTIC.  TAB 
Unannounced 
Justified  i  on, 


n 


Distribution/ 
Availability  Ccdes_ 
Avail  and/or 
Dist  !  Special 


n 


U 


i 


TR  6659 


LIST  OF  ILLUSTRATIONS 

Figure  Page 

1  Finite-Difference  Treatment  of  the  Horizontal  Interface  .  2-7 

2  Hierarchical  Structure  of  IFD  Model  .  3-2 

3  SVP  Tables .  3-5 

4  Summary  of  Bottom  Treatment  .  3-7/3-9 

5  Propagation  Loss  Versus  Range  for  Shallow  Water 

Propagation .  4-9 

6  Horizontal  Interface  Problem  .  4-10 

7  Solution  of  Horizontal  Interface  Problem  .  4-10 

8  Shallow-to-Deep  Water  Propagation  .  4-12 

9  Propagation  Loss  Versus  Range  for  Shallow-to-Deep  Water 

Propagation,  0.85  degree  Slope  .  4-13 

10  Propagation  Loss  Versus  Range  for  Shallow-to-Deep  Water 

Propagation,  8.5  degree  Slope  .  4-14 

11  Deep- to-Shal low  Water  Propagation  .  4-17 

12  Propagation  Loss  Versus  Range  for  Deep-to-Shallow 

Water  Propagation,  0.85  degree  Slope  .  4-18 

13  Propagation  Loss  Versus  Range  for  Deep-to-Shallow 

Water  Propagation,  8.5  degree  Slope  .  4-19 

14  Shallow-to-Deep  Water  Propagation,  Wedge-Shaped 

Region  With  a  Rigid  Sloping  Bottom  .  4-23 

15  Propagation  Loss  Versus  Range,  Wedge-Shaped 

Region  With  a  Rigid  Sloping  Bottom .  4-24 

Table 

1  Comparison  of  IFD  and  Exact  Solutions .  4-3 


11 


TR  6659 


IFD:  AN  IMPLICIT  FINITE-DIFFERENCE  COMPUTER 
MODEL  FOR  SOLVING  THE  PARABOLIC  EQUATION 


1.  INTRODUCTION 

The  parabolic  equation  (PE)  method  for  solving  a  class  of  range-dependent 
underwater  acoustic  wave  propagation  problems  and  the  split-step  Fourier 
algorithm  for  numerically  solving  the  parabolic  equation  were  introduced 
to  the  acoustics  community  by  Tappert  and  Hardin.1  Since  then,  several 
research  laboratories  have  implemented  the  split-step  method  into  various 
computer  programs  for  modeling  sound  propagation  in  the  ocean.  In  ocean 
environments  where  sound  interacts  weakly  with  the  bottom,  the  split-step 
method  is  fast  and  accurate,  but  in  cases  where  the  bottom  interaction  is 
strong,  the  method  is  less  accurate.  It  is  for  this  reason  that  the  search 
for  a  general  purpose  and  accurate  method  for  solving  the  parabolic  equation 
continues. 

Among  the  general  purpose,  accurate  numerical  methods  for  solving 
parabolic  wave  equations,  the  authors  have  found  two  classes  of  methods  to 
be  promising:  finite-differences  and  numerical  ordinary  differential 
equations.  In  this  report,  we  describe  a  finite-difference  approach  to 
solving  the  parabolic  equation. 

If,  as  the  PE  solution  is  marched  out  in  range,  the  boundary  information 
at  the  advanced  range  level  can  be  expressed  in  terms  of  the  known  values  at 
the  present  range  level,  then  implicit  finite-difference  (IFD)  methods  are 
more  desirable  than  explicit  finite-difference  methods  since  IFD  methods 
are  faster  and  unconditionally  stable.  For  these  reasons,  an  IFD  method 
for  solving  the  PE  was  selected  to  be  programmed  into  a  computer  model. 


The  IFD  model  presented  in  this  report  can  be  used  to  oredict  acoustic 
propagation  loss  in  both  range-dependent  and  range- independent  environments. 

An  important  feature  of  the  model  is  that  it  can  handle  arbitrary  surface 
boundary  conditions  and  an  irregular  bottom  with  arbitrary  bottom  boundary 
conditions.  In  the  event  that  the  bottom  boundary  conditions  cannot  be 
expressed  mathematically,  the  model  is  capable  of  introducing  artificial 
bottom  attenuation  to  make  the  problem  solvable.  Another  important  feature 
of  the  model  is  that  it  can  handle  horizontal  interfaces  of  layered  media. 

The  formulation  of  this  IFD  scheme  for  the  solution  of  the  parabolic 
wave  equation  is  discussed  in  the  next  section;  the  treatment  of  the  Neumann 
boundary  condition,  horizontal  interfaces,  and  attenuation  is  also  included. 
Section  3  presents  the  structure  of  the  IFD  model;  input  and  output  formats 
are  described  in  detail.  Section  4  is  devoted  to  test  problems;  input  run- 
streams  and  subroutines  are  included  in  the  event  the  reader  may  be  interested 
in  installing  the  IFD  model  in  his  computer. 

Program  listings  of  the  IFD  model  can  be  found  in  appendix  A.  Appendix  B 
contains  a  listing  of  a  program  which  may  be  used  to  plot  the  output  of  the 
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IFD  model;  this  program  is  annotated  and  self-explanatory.  All  programs  are 
written  in  FORTRAN  for  the  VAX-11/780  computer. 

It  is  requested  that  the  authors  be  notified  if  any  difficulties  with 
the  model  are  experienced.  User  contributions  that  will  enhance  the  model 
are  invited. 
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2.  SOLUTION  BACKGROUND 

The  theoretical  derivation  of  the  parabolic  equation  (PE)  and  the  IFD 
formula  for  solving  this  equation  are  summarized  below.  Details  concerning 
the  theoretical  development  of  the  PE  approximation  and  the  stability, 
consistency,  and  convergence  of  the  IFD  formula  can  be  found  in  references 
2  and  3.  Reference  3  also  contains  a  detailed  treatment  of  bottom  boundary 
conditions. 


2.1  THEORETICAL  DERIVATION 

We  begin  with  the  Helmholtz  equation  in  the  form 


v2P  +  k2  n2  P  =  0  ,  (2.1) 

where 

k  =  reference  wavenumber  =  -7- 
° 

n  =  n(r,z)  =  index  of  refraction  =  c^zy 
p  =  acoustic  pressure 
2 

v  =  Laplacian  operator 

cQ  =  reference  sound  speed 
c(r,z)  =  sound  speed 

u  =  2irf 

f  =  source  frequency. 

If  we  assume  a  geometry  that  is  cylindrically  symmetric,  equation  (2.1)  can 
be  written  as 


d~Z  +  r  +  ko  °2p  =  0  •  (2'2) 

dr  r  ar  0 

Using  the  parabolic  decomposition  technique,  p(r,z)  =  u(r,z)  v(r),  and 
neglecting  the  radial  dependence  field,  v(r),  we  obtain 

V  *  V  +  2,ko  ur  *  ko("2  *  0  “  ■  0  • 
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If  we  apply  farfield  and  paraxial  approximations,  the  second  derivative 
of  u  with  respect  to  r  may  be  dropped  from  the  above  equation,  which 
results  in  the  parabolic  equation  for  the  transmitted  field. 


iko<"  -  •>  .  t  u 

U - 2 - >'ra 

0 


(2.3) 


Equation  (2.3)  was  first  introduced  to  the  underwater  acoustics  community 
by  Tappert.4  Expressing  equation  (2.3)  in  a  general  form,  we  have 


u„  =  a(k  .r,z)u  +  b(k  ,r,z)u,_  , 

r  o  o  zz 


(2.4) 


where 


ik  (n^  -  1) 
a(kQ,r,z)  = - 2 - 


b(kQ,r,z)  =  2kQ 

Given  an  initial  field  and  the  appropriate  surface  and  bottom  boundary 
conditions,  an  initial  value  problem  of  equation  (2.4)  is  said  to  be  well 
posed  in  the  sense  of  Hadamard5  if  and  only  if  the  solution  exists,  is 
unique,  and  depends  continuously  on  the  initial  values. 

2.2.  IMPLICIT  FINITE-DIFFERENCE  FORMULA 

Using  the  Taylor  expansion,  we  obtain  a  two-level  scheme  for  equation 
(2.4): 


k_L 

=  e  3r  u(r,z) 


(2.5) 


If  we  let  z  =  mh,  r  =  nk,  and  u(r,z)  =  u(nk,mh)  =  u^  ,  where  the  depth 

and  range  increments  az  and  ar  are  denoted  by  h  and  k,  respectively,  then 
equation  (2.5)  can  be  written  as 


kJ_ 

un+1  =  e  3r  u" 
m  m 


(2.6) 
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(2.10) 

The  two  components  of  the  first  column  vector  on  the  right-hand  side  (RHS)  are 
the  two  boundary  points  at  the  advanced  range  level.  The  two  components  of 
the  last  column  vector  on  the  RHS  are  the  two  boundary  points  at  the  initial 
range  level.  This  scheme  is  known  as  the  Crank-Nicolson6  scheme  for  variable 
coeffi cents. 

In  a  region  where  the  reference  sound  speed  is  constant,  bm  is  also 
constant.  Where  bm  is  constant,  equation  (2.10)  may  be  simplified  such  that 
the  matrix  form  may  be  expressed  as 

x  u"*1  -  ♦ » un  ♦  uj  . 


2-4 


where 


X  is  a  square  matrix  whose  elements  are 


n+1  _  Yn+1 

i,i+l  "  Ai+l,i  ‘  " 


0  elsewhere. 

un+1  is  a  vector  whose  components  are  the  unknowns  u"+1  at  the  advanced 
level.  1 

u"+1  is  a  vector  having  zero  components  everywhere  except  for  the  first 
and  last  components,  which  are  the  surface  and  bottom  boundary  conditions, 
respectively,  at  the  advanced  range  level. 


Y  is  a  square  matrix  whose  elements  are 


m 


y n  _  Yn 
i ,i+l  ”  Ti+l,i 


=  1 


0  elsewhere. 

un  is  a  vector  whose  components  are  the  known  values  u"  at  the  present 
range  level . 

u*  is  a  vector  having  zero  components  everywhere  except  for  the  first 
and  the  last  components,  which  are  the  surface  and  bottom  boundary 
conditions,  respectively,  at  the  present  range  level. 

The  index  i  ranges  from  1  through  m. 


2.3  BOUNDARY  TREATMENT 

An  optional  homogeneous  Neumann  boundary  condition  that  has  been  pro¬ 
grammed  into  the  package  is  described  below.  However,  when  other  types  of 
boundary  conditions  arise,7  the  user  may  program  an  optional  subroutine  that 
will  allow  the  program  to  accept  the  existing  boundary  conditions. 

Consider  p^  =  0,  which  gives 

P2  cos  0  -  pr  sin  0  =  0  , 


L 
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where  0  is  the  angle  of  the  bottom  with  respect  to  the  surface.  Since  u^ 
satisfies  our  parabolic  wave  equation,  then  the  associated  boundary  condition 
for  equation  (2.4)  can  be  described  by  the  second-order  ordinary  differential 
equation 


cot  0 


zz 


u2  + 


»"o1)(|y> 

hTtt^H 


Fu  s  o 


(2.11) 


Expressing  equation  (2.11) 


by  a  second  order  finite-difference,  we  have 


>  -  *  “PH*  ♦  ‘ -s 


c0-  -  -2  +  h2(1k„  +  a)  £)«„  +  Vl-° 


Solving  the  above  equation  for  um+1  in  terms  of  um  and  um_^,  we  find  that 


um+l  Q  um  "  Q  um-l 


(2.12) 


where 


P  .  h  SJti  .  2  ♦  h2(1k  *  a)  i 


n  _  i  u  COt  0 

Q  _  !  .  h  — g—  . 

The  coefficients  of  um  and  um-i  are  combined  with  the  appropriate  Y  and  X 
on  both  sides  of  equation  (2.10).  This  mathematical  manipulation  gives  us 
the  advantage  of  using  an  unconditionally  stable  IFD  formula.2 


2.4  INTERFACE  TREATMENT 

Included  in  the  IFD  model  is  the  finite-difference  treatment  of  the 
horizonta1  interface,  as  shown  in  figure  1.  The  densities,  pi  and  P2,  are 
considered  to  be  constant  throughout  the  present  treatment,  to  be  consistent 
with  the  finite-difference  formulation,  az  is  denoted  by  h,  and  ar  is  denoted 
by  k.  Z3  is  the  depth  of  the  interface  boundary.  The  subscripts  1  and  2 
denote  medium  1  and  medium  2,  respectively.  As  shown  in  reference  8,  the 
horizontal  interface  parabolic  equation  is 


(2.13) 
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Figure  1.  Finite-Difference  Treatment  of  the  Horizontal  Interface 


Using  the  same  implicit  finite-difference  scheme  to  solve  equation  (2.13),  we 
have 


where 


i  -  j  k  p"+1  Qr1  -  k  p"+1 
2mm  m 


xn+1)  u 
ZZ  /  I 


n+1 

m 


1  +  i  k  P;  Q"  +  k  Pm 
2mm  m 


n  n  \  n 
m  Tzz)  um  * 


P-lf  + 

D,  p0  0, 


v-1 


'1  y2  2> 


i  '1  P1  a2 

Q  =  (  r-  +  —  r- 
\bl  p2  z  / 


(2.14) 


(2.15) 


(2.16) 


Tzz  *  u  =  ^  Jl  (Vl  '  um)  ’  (um  •  um-l) 


(2.17) 
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In  the  original  IFD  mathematical  formulation,2  the  m-th  equation  reads 


1  Jc  hn+l  .  1  .  n+1  k  .  n+1  1  k  .  n+ll 

2~Z  bn,  •  1  *  2  ka«  ■  T-2  bm  , 


n+1 

m 

n+1 

m+1 


j 


1  k 


n 


*  ^  u 

2  7?  ^m’ 

^  h 


It  ±ka"-i  b". 
2  m  ^2  m 


1  k  Kn 


r  n  ^ 
Vi 

n 

um 

n 

vum+l 


(2.18) 


Writing  equation  (2.14)  in  matrix  form,  we  have 


\  Pm+1*  1 
h2  m 


j  k  p"+1  q"+1  +  4  Pn+1 

2  m  m  ^4 


m 


i*r  • 

V  °2/ 


4  C1  -i 

h2  "  p-2i 


k  0n 
m‘ 


P" 


1  +  7  k  Pm  Qm  ’  4  Pm  l1  +  ~ ]  • 

2  mm  |^4  mV  P2  / 


4  Pn  !i 

h2  m  p2 


c? 

C1 


n+i 

'm+l 


r«tr 

n 

‘m 
n 

lVu 


u 


(2.19) 


Since  P  is  range  independent,  equation  (2.19)  may  be  simplified  such 
that  we  have 
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r..n+l"\ 

Jm-1 


■1,  y-  (P)"1  -  |  h2  Qn+1  + 

m 


i,  nT  (p )■*  ♦  *  h‘  q; 


u"*1 
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u!!+1 

.  nH-l 


"cr 


m 

Vi 

V.  J 


(2.20) 


Incorporating  the  horizontal  interface  treatment  into  the  present  IFD  model 
requires  replacing  the  row  vector  elements  of  both  sides  of  equation  (2.18) 
by  the  corresponding  row  vector  elements  of  equation  (2.20). 

Rewriting  equation  (2.10),  we  have 
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where 


r  -  r  w1  -  t  »2  c1  + 


:  -  V  (P)-1  -  i  Qj 


The  IFD  model  presented  in  this  report  solves  equation  (2.21). 


2.5  ATTENUATION 

Attenuation  is  applied,  as  suggested  by  Jensen  and  Krol,9  by  inserting 
the  units  of  loss  in  the  imaginary  part  of  the  index  of  refraction  as  shown 
below: 

2  _/co  f  +  ,  (Co  \  6 

n  '[cTj  ^c7  1  27.  287527 

where 


S  is  the  attenuation  in  dB/wavelength 
cQ  is  the  reference  sound  speed  in  m/s 
c^  is  the  speed  of  sound  in  m/s  at  depth  i. 
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3.  COMPUTER  IMPLEMENTATION 

The  IFD  model  that  implements  the  implicit  finite  difference  formula, 
equation  (2.21),  has  been  written  in  FORTRAN  using  single  precision,  complex 
arithmetic  and  has  been  installed  on  a  VAX-11/780  digital  computer.  A  listing 
of  the  model  can  be  found  in  appendix  A.  The  program  listing  is  heavily 
"commented"  in  the  event  the  reader  wishes  to  trace  through  the  program  logic. 


3.1  DESIRABLE  FEATURES 

The  desirable  features  to  be  noted  in  the  IFD  model  are  classified  below. 

A.  General i ty 

The  model  can  solve  the  parabolic  equation  in  the  general  form 

ur  =  a(V  r’z^  u  +  b0<o’r*2)  uzz  ' 

The  following  environments  can  be  handled: 

•  Range-dependent/ range- i ndependent 

•  Shallow  water/deep  water 

•  Shallow-to-deep  water,  deep-to-shallow  water,  or  the  combination. 

B.  Portabi 1 i ty 

All  programs  are  written  in  FORTRAN  and  can  be  installed  In  most 
other  computers  without  much  difficulty. 

C.  Flexibility  and  Reliability 

The  user  may,  at  his  option,  include  his  own  versions  of  the  sub¬ 
routines  that  input  environmental  parameters. 

In  the  event  the  bottom  boundary  information  is  unavailable,  our  problem, 
equation  (2.4),  becomes  111  posed.  Theoretically,  this  problem  is  not 
uniquely  solvable  and  the  computation  should  not  be  performed.  However, 
by  applying  the  artificial  bottom  technique,  the  bottom  may  be  extended  deep 
enough  such  that  a  zero  boundary  condition  results,  whereupon  the  problem 
becomes  artificially  well  posed  and  a  unique  solution  exists. 

D.  Accuracy 

The  Initial  local  truncation  error2  of  the  IFD  model  is 
E( I]  -  0(k3  ♦  kh2)  . 
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Accuracy  comes  automatically  from  theory  and  can  be  controlled  by  a  PC 
(predictor-corrector)  technique  at  the  expense  of  computation  overhead 
costs.  The  present  version  of  the  IFD  model  does  not  include  PC  techniques. 

E.  Special  Features 

Special  features  include 

•Automatic  handling  of  an  irregular  bottom 

•Acceptance  of  arbitrary  bottom  boundary  conditions 

•Automatic  handling  of  a  homogeneous  Neumann  bottom  boundary 
condition 

•Automatic  handling  of  multiple  horizontal  interfaces  and  layers 
•  Capability  to  introduce  an  artificial  absorbing  layer. 


3.2  STRUCTURE  OF  PROGRAM 

The  IFD  computer  model  consists  of  a  main  program  and  10  subroutines. 
The  hierarchical  structure  of  the  model  is  as  shown  in  figure  2. 


Figure  2.  Hierarchical  Structure  of  IFD  Model 
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Those  subroutines  marked  with  an  asterisk  (*)  are  prepared  by  the  user; 
the  subroutine  marked  with  a  dagger  (f)  may  be  modified  by  the  user.  Input 
parameters  that  result  in  control  being  transferred  to  user-supplied  sub¬ 
routines  are  as  shown  in  the  figure.  For  example,  if  input  parameter  ISF*1, 
then  control  is  transferred  to  UFIELD  rather  than  SFIELD. 

Linking  the  program  is  performed  as  follows:  LINK  IFD,  CRNK,  TRID, 

DIAG,  HNKL,  SVP,  SFIELD,  USVP,  UFIELD,  BCON,  SCON.  A  brief  description  of 
each  subroutine  follows. 


3.2.1  Main  Program  IFD 

IFD  is  the  main  program  of  the  model  and  controls  execution  of  the 
various  subroutines  which  make  up  the  model.  Initially,  IFD  reads  selected 
input  parameters  and  performs  initialization  of  certain  variables.  IFD  then 
calls  on  either  subroutine  SFIELD  or  UFIELD  to  generate  the  starting  field 
that  is  to  be  marched  out  in  range.  Selected  problem  parameters  are  then 
printed  and,  if  requested,  written  in  an  output  file  for  subsequent  use. 

IFD  then  calls  on  subroutine  DIAG  to  compute  the  main  diagonals  of  the 
matrices  that  represent  the  system  of  equations  at  the  present  and  advanced 
ranges . 

After  these  preliminary  procedures  have  been  accomplished,  IFD  enters 
a  main  loop  and  continues  to  cycle  in  the  loop  until  the  solution  has  been 
marched  out  to  maximum  range  as  requested. 

While  in  the  main  loop,  as  the  solution  is  marched  out  in  range,  IFD 
determines  whether  or  not  to  update  the  sound  speed  profile  and/or  bottom 
depths.  If  an  update  is  performed,  IFD  calls  on  subroutine  DIAG  to  recompute 
the  main  diagonals  in  the  matrices.  Whether  or  not  the  diagonals  have  been 
updated,  IFD  calls  on  subroutine  CRNK  to  advance  the  solution  one  range  step, 
aR,  forward.  The  solution  returned  by  CRNK  is  then  printed  and  written  in 
an  output  file  as  requested  by  the  user.  When  the  solution  range  has  reached 
the  maximum  range,  the  program  is  terminated.  If  the  solution  range  has  not 
reached  the  maximum  range,  IFD  returns  to  the  top  of  the  main  loop  and  repeats 
the  above  procedures. 


3.2.2  Subroutine  SFIELD 


If  input  parameter  ISF  *  0,  main  program  IFD  calls  on  subroutine  SFIELD 
to  generate  a  Gaussian  starting  field  at  zero  range.  The  Gaussian  starting 
field  defined  by  Brock10  is 


where 


U(I)  *  CMPLX  (PR,  0.0)  , 


PR  =  GA 


e 


-e 
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ZM 

ZS 

GW 

FK 

GA 


depth  of  point  in  mesh  in  meters  (I*az) 

source  depth  in  meters 
2 


Gaussian  width  =■ 


FK 


reference  wavenumber 

1  /  ?  Vs 

Gaussian  amplitude  =  ^  l-^l 


3.2.3  Subroutine  UFIELD 


If  input  parameter  ISF  =  1,  main  program  IFD  calls  on  user-written  sub¬ 
routine  UFIELD  to  generate  the  starting  field.  UFIELD  must  load  N  values  of 
the  complex  starting  field  in  array  U. 

If  ISF  =  0,  UFIELD  is  not  called  and  may  be  a  dummy  subroutine. 

If  the  user-supplied  starting  field  is  an  elliptic  solution,  then  the 
starting  field  must  be  divided  by  the  Hankel  function  and  the  solution  field 
multiplied  by  the  Hankel  function  before  computing  transmission  loss.  TRe 
above  procedure  can  be  accomplished  automatically  by  setting  input  parameter 
IHNK  =  1. 


3.2.4  Subroutine  SVP 


When  the  range  of  the  solution  is  equal  to  the  range  of  the  next  sound 
velocity  (speed)  profile  (input  parameter  RSVP),  a  new  sound  speed  profile 
is  inputted.  If  input  parameter  KSVP  a  0,  then  subroutine  SVP  is  called 
upon  to  read  the  next  sound  speed  profile  from  the  input  runstream.  SVP  then 
loads  the  six  data  tables  shown  in  figure  3. 

Table  ZLYR  contains  the  maximum  depth  of  each  layer  as  measured  from 
the  surface.  Tables  RHO  and  BETA  contain  the  density  and  attenuation, 
respectively,  in  each  layer.  Variable  NLYR  is  the  number  of  layers  inputted. 
The  sound  speed  profile  for  each  layer  is  stored  in  tables  ZSVP  and  CSVP. 

NSVP  is  the  total  number  of  points  stored  in  these  tables.  If  an  error  is 
detected  while  inputting  the  profiles,  NSVP  is  set  to  0.  Table  IXSVP  contains 
indexes  that  point  to  the  last  sound  speed  profile  value  in  each  layer. 

At  the  present  stage  of  development,  linear  interpolation  of  sound  speed 
values  is  performed  only  in  depth.  Changes  in  the  profiles  in  range  are 
abrupt,  with  no  interpolation  performed.  Interpolation  of  the  sound  speed 
values  is  performed  in  subroutine  DIAG. 
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ZLYR 


RHO 


BETA 


MAXIMUM  DEPTH 

OF  LAYER  1 

DENSITY  IN 
LAYER  1 

ATTENUATION 

IN  LAYER  1 

MAXIMUM  DEPTH 

OF  LAYER  2 

DENSITY  IN 
LAYER  2 

ATTENUATION 

IN  LAYER  2 

1 

MAXIMUM  DEPTH 

OF  LAYER  NLYR 

DENSITY  IN 
LAYER  NLYR 

ATTENUATION 

M  LAYER  NLYR 

*  EXTENDED  IF  USER  REQUESTS  AN  ARTIFICIAL  ABSORBING  LAYER. 


Figure  3.  SVP  Tables 
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3.2.5  Subroutine  USVP 

If  input  parameter  KSVP  *  0,  then  subroutine  USVP  is  called  to  supply 
an  updated  SVP  at  each  step  in  range.  Subroutine  USVP  must  be  prepared  by 
the  user,  who  must  load  the  six  tables  described  in  the  previous  section. 
Variables  NLYR  and  NSVP  must  also  be  loaded  by  the  user. 

Variable  KSVP  may  be  used  in  a  computed  GOTO  statement  to  transfer 
control  within  user  program  USVP.  When  the  user  no  longer  needs  USVP, 

KSVP  must  be  set  to  zero  within  USVP.  The  last  profile  entered  will  be  used 
until  the  solution  range  is  equal  to  the  next  RSVP.  If  KSVP  is  not  set  to 
zero,  then  USVP  will  be  called  until  the  solution  range  is  equal  to  RSVP, 
the  range  of  the  next  profile.  With  this  option,  the  user  can  generate  a 
new  profile  at  each  range  step.  Sound  speed  profile  values  interpolated  in 
range  may  be  entered  by  the  user  in  this  manner.  When  the  next  solution 
range  is  equal  to  the  next  RSVP,  either  SVP  or  USVP,  depending  on  the  next 
KSVP,  is  called  to  input  the  next  profile. 


3.2.6  Subroutine  DIAG 

Subroutine  DIAG  computes  the  range-dependent  and  depth-dependent  main 
diagonals  of  the  matrices  that  represent  the  system  of  equations  at  the 
present  and  advanced  ranges,  as  shown  in  equation  (2.21). 

Prior  to  computing  the  diagonals,  DIAG  determines  the  values  of  sound 
speed,  density,  and  attenuation  to  be  used  at  each  depth  represented  by 
the  corresponding  row  of  each  matrix.  Linear  interpolation  of  sound  speed 
in  depth  is  performed  as  required. 


3.2.7  Subroutine  CRNK 

Subroutine  CRNK  computes  the  right-hand  side  of  the  system  of  equations, 
D(i);  determines  bottom  type;  sets  up  bottom  conditions  at  the  present  and 
advanced  ranges;  and  then  calls  on  subroutine  TRID  to  solve  the  tridiagonal 
system  of  equations.  If  the  user  is  supplying  surface  conditions,  then  CRNK 
calls  on  SCON  to  provide  SURY  and  SURX,  which  are  added  to  the  computation  of 
D(l).  SURY  and  SURX  are  the  values  of  the  field  at  the  surface  at  the  present 
and  advanced  ranges,  respectively.  If  the  user  is  supplying  bottom  conditions, 
then  CRNK  calls  on  user-written  subroutine  BCON  for  these  conditions.  A  sum¬ 
mary  of  the  treatment  of  the  bottom  follows  in  figure  4. 

Variables  U ( 1 )  through  U(N)  represent  the  known  values  of  the  field  at 
the  present  range.  Variables  U(l)'  through  U(N)'  represent  the  values  of  the 
field  at  the  advanced  range  and  are  to  be  determined.  D(N)  is  the  right-hand 
side  of  row  N  in  the  system  of  equations.  D(l)  through  D(N-l)  are  computed 
prior  to  performing  the  operations  summarized  on  the  following  pages.  BOTY 
and  BOTX  are  the  values  of  the  field  at  the  bottom  at  the  present  and  advanced 
ranges,  respectively.  ITYPEB  is  an  input  parameter  for  selecting  the  type  of 
bottom  treatment  to  be  performed. 
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ITYPEB  -  0,  RIGID  BOTTOM,  FLAT,  9  *  0.0C 
AR 

SURY  ( - * - \  ^SURX 

U(1)' 


U(N) 


BOTY 


U(N)' 


BOTX 


} 


AZ 


SUMMARY 
SET  BOTY  -  U(N) 

COMPUTE  D(N) 

ASSUME  BOTX  -  U(N)' 
SOLVE  FOR  U(1)'  THRU  U(N)' 


ITYPEB  =  0,  RIGID  BOTTOM,  SLOPES  UP,  9 <0.0° 
SURY  SURX 


•  •  U(N —  1 )' 


SUMMARY 
SET  BOTY  -  U(N) 

SET  N  »  N  —  1 
COMPUTE  D(N) 

ASSUME  BOTX  =*  f  (U(N)',  U(N-1)',  V  V  . VK) 

SOLVE  FOR  U(1 )'  THRU  U(N)' 


ITYPEB  -  0,  RIGID  BOTTOM,  SLOPES  DOWN,  9  >0.0° 


SURY 


SURX 


U(N  — 1)  •  U(N—  1 )'  • 


SUMMARY 

ASSUME  BOTY  -f(U(N),  U(N-1),  Vr  Vg . ,VK) 

COMPUTE  D(N) 

ASSUME  BOTX  «f(U(N)\  U(N-1)',  V  V  . VK) 

SOLVE  FOR  U(1 )'  THRU  U(N)' 

SET  N  -  N  +  1 

U(N)'  -nU(N-l)’,  U(N  — 2)',  V  V  . VK) 


BOTX,  NEW  U(N)' 

Summary  of  Bottom  Treatment 
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ITYPEB  -  1,  USER  BOTTOM,  FLAT,  O.O8 


SURY 


SURX 


SUMMARY 

BOON  SETS  BOTY  -  U(N) 
BOON  SUPPLIES  BOTX 
N  -  N  -  1 
COMPUTE  D(N) 

SOLVE  FOR  U(1 )'  THRU  U(N)' 

N  «  N  +  1 

SET  U(N)'  -  BOTX 


ITYPEB  =  1 ,  USER  BOTTOM,  SLOPES  UP,  0<O.O® 


SURY  SURX 


■  • 


SUMMARY 

BOON  SETS  BOTY  -  U(N) 
BOON  SUPPLIES  BOTX 
N  -  N  -  1  (DELETE  A  POINT) 
COMPUTE  D(N) 

SOLVE  FOR  U(1 )'  THRU  U(N) 


ITYPEB  «  1 ,  USER  BOTTOM,  SLOPES  DOWN,  9  >0.0° 


SURY  SURX 


SUMMARY 

BCON  SUPPLIES  BOTY  AND  BOTX 
COMPUTE  D(N) 

SOLVE  FOR  U(1 )'  THRU  U(N)' 

N  -  N  +  1  (ADD  A  POINT) 

SET  U(N)'  -  BOTX 


Figure  4b.  Summary  of  Bottom  Treatment 


3-8 


ITYPEB  -  2,  ARTIFICIAL  ABSORBING  LAYER,  ~l_AT  9  -  0.0° 
SEE  ITYPEB  =*  3 


ITYPEB  -  2,  ARTIFICIAL  ABSORBING  LAYER.  SLOPES  UP,  0<O.O# 


SURY  SURX 
— » — — 


U(N) 


ARTIFICIAL 

LAYER 


SUMMARY 
U(N)  -  0.0 
Set  U(N— 1)'  -  0.0 
Set  U(N)‘  -  0.0 

SOLVE  FOR  U(1)'  THRU  U(N-2)# 
N  -  N  -  1  (DELETE  A  POINT) 


ITYPEB  -  2,  ARTIFICIAL  ABSORBING  LAYER,  SLOPES  DOWN,  0>O.O° 


SURY  SURX 


SUMMARY 

*  *  U(N)  -  BOTY  -  BOTX  -  0.0 

COMPUTE  D(N) 

•  •  SOLVE  FOR  U(1 )'  THRU  U(N)' 

N-N+1  (ADD  A  POINT) 

Set  U(N)'  -  0.0 


ARTIFICIAL 

LAYER 


ITYPEB  -  3.  ARTIFICIAL  ABSORBING  LAYER,  FLAT,  9  -  0.0° 


SURY  SURX 


SUMMARY 
U(N)  -  0.0 
SET  U(N)'  -  0.0 

SOLVE  FOR  U(1)'  THRU  U(N-1)' 


ARTIFICIAL 

LAYER 


Figure  4c.  Summary  of  Bottom  Treatment 
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3.2.8  Subroutine  BCON 


If  input  parameter  ITYPEB  =  1,  then  subroutine  CRNK  calls  on  user- 
written  subroutine  BCON  to  supply  BOTY  and  BOTX,  the  values  of  the  field 
on  the  bottom  at  the  present  and  advanced  ranges,  respectively.  The  treat¬ 
ment  of  the  user  bottom  was  described  earlier  in  the  description  of  sub¬ 
routine  CRNK. 


3.2.9  Subroutine  SCON 


If  the  user  wishes  to  supply  values  for  SURY  and  SURX,  the  values  of 
the  field  at  the  surface  at  the  present  and  advanced  ranges,  respectively, 
he  may  do  so  by  rewriting  subroutine  SCON.  At  present,  subroutine  SCON 
sets  SURY  and  SURX  to  0.0. 


3.2.10  Subroutine  TRIO 

Subroutine  TRIO,  a  specialized  version  of  subroutine  TRIDAG  presented 
in  reference  11,  solves  a  system  of  N  linear  simultaneous  equations  having 
a  tridiagonal  coefficient  matrix.  As  shown  in  equation  (2.21),  the  lower 
diagonal  coefficient  at  the  advanced  range  in  rows  2  through  N  is  -1.  The 
upper  diagonal  coefficient  in  rows  1  through  N-l  is  the  negative  of  the 
ratio  of  the  densities  of  the  media  above  and  below  the  receiver  depth 
represented  by  each  row.  If  a  receiver  lies  entirely  within  the  same  medium, 
then  the  upper  diagonal  coefficient  for  that  row  is  also  -1.  The  main 
diagonal  is  computed  in  routine  DIAG. 


3.2.11  Complex  Function  HNKL 

HNKL  computes  the  Hankel  function.  The  algorithm  for  computing  the 
Hankel  function  is  described  in  reference  12. 


3.3  USER'S  GUIDE 

The  next  two  subsections  describe  input  and  output  formats  in  detail. 
Test  examples  showing  sample  input  runstreams  and  user-written  subroutines 
can  be  found  later  in  this  report. 

Since  plotting  facilities  vary  from  one  activity  to  another,  it  was 
decided  that  the  plotting  program  should  not  be  included  in  the  IFD  model. 
A  separate  plot  program  has  been  included  in  appendix  B  for  reference 
purposes. 


3.3.1  Input  Format 

Prior  to  executing  the  IFD  model,  input  card  images  containing  problem 
parameters  must  be  stored  in  file  IFD. IN.  File  IFD. IN  is  assigned  to  FORTRAN 
unit  number  NIU  in  the  main  program.  If  the  user  prefers  to  input  problem  param¬ 
eters  on  cards,  then  parameter  NIU  should  be  equated  to  the  card  reader  unit 
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number,  and  the  statement  which  assigns  file  IFD.IN  should  be  removed 
from  the  main  program.  In  either  case,  the  input  runstream  is  prepared 
in  free  format  as  follows: 

CARD  CONTENTS 

1  FRQ,  ZS,  CO,  ISF,  RA,  ZA,  N,  IHNK,  I TYPES,  ITYPES 

where 


FRQ  =  frequency  (Hz) 

ZS  =  source  depth  (m) 

CO  =  reference  sound  speed  (m/s) 

If  CO  =  0.0,  CO  is  set  to  the  average  sound  speed  in 
the  first  layer. 

ISF  =  starting  field  flag. 

0  =  Gaussian  starting  field  is  generated. 

1  =  user  prepares  starting  field.  See  subroutine  UFIELD. 

If  ISF  =  0,  RA  is  set  to  zero. 

RA  =  horizontal  range  from  source  to  starting  field  (m). 

If  ISF  =  0,  RA  is  set  to  0.0. 

ZA  *  depth  of  starting  field  at  range  RA  (m).  If  ZA  *  0.0, 

ZA  is  set  to  the  maximum  depth  of  the  bottom-most 
layer  in  the  first  profile. 

If  ITYPEB  =  2  or  3  and  ZA  *  0.0,  ZA  is  set  to  (4/3)* 
maximum  depth  of  the  bottom-most  layer.  If  ITYPEB 
=  2  or  3  and  ZA  +  0.0,  the  artificial  bottom  layer  is 
extended  to  ZA  meters  provided  that  ZA  is  greater 
than  or  equal  to  the  maximum  depth  of  the  bottom-most 
layer  in  the  first  profile. 

N  *  Number  of  equi spaced  receivers  in  the  starting  field. 

If  N  =  0,  N  is  set  so  that  the  receiver  depth  increment 
is  less  than  or  equal  to  1/4  wavelength.  If  N  is  greater 
than  MXN,  N  is  set  to  MXN.  See  parameter  MXN. 

IHNK  *  Hankel  function  flag.  If  IHNK  *  0,  don't  use  Hankel 
function.  If  IHNK  =  1,  divide  the  starting  field  by 
the  Hankel  function,  then  multiply  the  solution  field 
by  the  Hankel  function  before  computing  propagation 
loss.  If  starting  field  is  Gaussian,  IHNK  should  be 
set  to  0.  If  starting  field  is  ell iptic,  IHNK  should 
be  set  to  1. 
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CARD 


2 

where 
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CONTENTS 

ITYPEB  =  type  of  bottom 

s  0,  homogeneous  Neumann  boundary  condition: 
program  supplies  bottom  condition 

3  1,  user  supplies  bottom  condition.  See  subroutine 
BCON. 

=  2,  artificial  absorbing  layer  introduced:  bottom 
of  layer  follows  contour  of  water-bottom  interface. 

=  3,  artificial  absorbing  layer  introduced:  bottom 
of  layer  kept  flat. 

ITYPES  3  type  of  surface 

3  0,  pressure  release:  SCON  sets  SURY  and  SURX  3  0.0. 

/  0,  user  inserts  code  in  SCON  to  compute  SURY  and 
SURX. 

RMAX,  DR,  WDR,  WDZ,  PDR,  PDZ,  ISFLD,  ISVP,  IBOT 

RMAX  3  maximum  range  of  solution  (m). 

DR  3  range  step  for  marching  solution  (m). 

If  DR  3  0,  DR  is  set  to  1/2  wavelength. 

WDR  3  range  step  (rounded  to  nearest  DR)  at  which  solution 
is  written  on  disk  (m). 

WDZ  =  depth  increment  (rounded  to  nearest  DZ)  at  which 
solution  is  written  on  disk  (m). 

PDR  *  range  step  (rounded  to  nearest  DR)  at  which  solution 
is  printed  (m). 

PDZ  3  depth  increment  (rounded  to  nearest  DZ)  at  which 
solution  is  printed  (m). 

ISFLD  3  0,  don't  print  starting  field. 

3  1,  print  starting  field. 

ISVP  *  0,  don't  print  sound  speed  profile. 

3  1,  print  sound  speed  profile. 

IBOT  3  0,  don't  print  bottom  depths. 

1,  print  bottom  depths. 
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CARD 

3 

4 

5 

N 

N+l 

N+2 

N+3 

N+4 

N+5 

N+6 

N+7 

N+M 


CONTENTS 
Rl,  Z1 
R2,  Z2 
R3,  Z3 


-1,  -1 


BOTTOM  PROFILE 

range  and  depth  of  water  (m). 
maximum  number  of  depths  =  100 
{See  parameter  MXTRK.) 


Required.  Marks  end  of  bottom  profile. 


RSVP 

KSVP 


NLYR 


ZLYR(I),  RHO(I),  BETA(I) 
ZSVP(l),  CSVP(l) 

ZSVP(2),  CSVP(2) 

ZSVP(J),  CSVP(J) 


REPEAT 
FOR  EACH 
LAYER. 
1=1,  NLYR 


REPEAT 
FOR  EACH 
PROFILE 


where 

RSVP  =  range  of  SVP  (m). 

KSVP  =  SVP  flag. 

3  0,  profile  is  in  runstream. 

t  0,  profile  (cards  N+4  through  N+M)  is  supplied  by 
user-written  subroutine  USVP.  KSVP  may  be  used 
in  computed  GOTO  statement  to  transfer  control 
in  user  subroutine  USVP. 

NLYR  3  number  of  layers.  If  ITYPEB  *  2  or  3,  the  program 

inserts  an  artificial  absorbing  layer  and  then  increments 
NLYR  by  1.  Maximum  NLYR  3  100  (see  parameter  MXLYR). 


3-13 


TR  6659 


ZLYR(I)  *  maximum  depth  of  layer  I  in  profile  (m). 

RHO(I)  =  density  In  layer  I  ig/cm^). 

BETA(I)  *  attenuation  in  layer  I  (dB/wavelength) 

ZSVP  -  SVP  depth  (»). 

Maximum  number  of  sound  speed  profile  values  3  100 
(see  parameter  MXSVP). 

ZSVP(l)  *  depth  to  top  of  layer  I  (m). 

ZSVP(O)  =  depth  to  bottom  of  layer  I  (m). 

CSVP  3  SVP  speed  (m/s). 

Maximum  number  of  sound  speed  profile  values  *  100 
(see  parameter  MXSVP). 

CSVP(l)  =  speed  of  sound  at  top  of  layer  I  (m/s) 

CSVP(J)  =  speed  of  sound  at  bottom  of  layer  I  (m/s). 


3.3.2  Output  Format 

Output  from  the  IFD  model  is  written  on  disk  in  file  IFD.OUT.  File 
IFD.OUT  is  assigned  to  FORTRAN  unit  number  NOU  in  the  main  program.  The 
data  written  in  IFD.OUT  are  unformatted  and  are  written  with  FORTRAN  WRITE 
statements  as  follows: 


WRITE(NOU)FRQ,ZS,CO,ISF,RA,ZA,N,IHNK,ITYPEB,ITYPES,RMAX,DR,WDR,DZ,NLYR,ZLYR,RHO,BETA 
WRITE(NOU)NZ,RA,WDZ,(U( I ) ,I=IWZ,N,IWZ) 


The  first  WRITE  statement  is  executed  only  once  and  writes  the  value 
of  each  of  the  following  parameters  at  the  start  of  the  problem. 

FRQ  *  frequency  (Hz) 

ZS  3  source  depth  (m) 

CO  3  reference  sound  speed  (m/s) 

ISF  3  starting  field  flag 

3  0,  Gaussian 

3  1,  user 
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RA  3  horizontal  range  from  source  to  starting  field  (m). 

ZA  =  depth  of  starting  field  (m). 

N  *  number  of  equispaced  receivers  in  the  starting  field. 

IHNK  *  Hankel  function  flag. 

3  0,  Hankel  function  not  used. 

=  1,  starting  field  was  divided  by  Hankel  function. 

ITYPEB  3  type  of  bottom. 

3  0,  rigid. 

=  1,  user  bottom. 

3  2,  artificial  absorbing  layer— follows  bottom  contour. 

=  3,  artificial  absorbing  layer-bottom  kept  flat. 

I TYPES  =  type  of  surface 

*  0,  pressure  release. 
t  0,  user  supplied. 

RMAX  =  maximum  range  of  solution  (m). 

DR  =  range  step  for  marching  solution  (m). 

WDR  =  range  step  at  which  solution  is  written  on  disk  (m). 

DZ  =  depth  increment  of  receivers  (m). 

NLYR  3  number  of  layers. 

ZLYR  3  array  containing  depth  of  each  layer  (m). 

3 

RHO  3  array  containing  density  in  each  layer  (g/cm  ). 

BETA  3  array  containing  attenuation  in  each  layer  (dB/ wave length) . 
The  second  WRITE  statement  is  executed  at  each  write-range  increment, 

WDR.  The  data  written  are  as  follows: 

NZ  =  number  of  equispaced  receivers  in  the  solution  field. 

RA  3  horizontal  range  from  source  to  solution  (m). 
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WDZ  =  depth  increment  at  which  solution  is  written  on  disk  (m) 
U  *  array  that  contains  the  complex  field  at  **ange  RA. 

If  IHNK*1,  then  the  contents  of  II  must  be  If’-l  ed 

by  the  Hankel  function  before  computing  propagati-  loss 

IWZ  *  Index  increment  of  receiver  solutions  to  be  wr  n 

disk. 

The  following  READ  statement  may  be  used  to  read  the  solution  field: 

READ (unit)  NZ,RA,WDZ, (U( I) ,I=1,NZ). 
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4.  TEST  PROBLEMS 

The  first  test  problem  is  an  exact  solution  that  has  no  physical 
significance  other  than  to  demonstrate  the  accuracy  and  flexibility  of  the 
IFD  model.  This  problem  requires  the  user  to  write  the  subroutines  that 
supply  the  starting  field,  sound  speed  profiles,  bottom  condition,  and 
surface  condition. 

Several  other  test  problems  that  demonstrate  the  capability  of  the  IFD 
model  are  also  included.  One  problem  deals  with  a  range- independent  environ¬ 
ment;  the  others  deal  with  range-dependent  environments.  Input  runstreams 
and  user-written  subroutines  that  produced  the  IFD  solutions  to  these  problems 
are  also  included. 


4.1  EXACT  SOLUTION 

As  a  test  for  accuracy,  the  IFD  model  was  used  to  solve  Burger's13 
kinetics-diffusion  parabolic  partial  differential  equation  for  which  an 
exact  solution  is  known.  Burger's  equation  is 


ut  =  vuxx  -uV°-x-1,t-°  ’ 


where  the  subscripts  denote  partial  derivatives.  An  exact  solution  for 
Burger's  equation  is 


u(x,t)  = 


[■ 


1  +  exp  (x/2v)  -  (t/4v) 


-1 


Substituting  r  for  t  and  z  for  x  and  rewriting  Burger's  equation  in  the 
general  form  of  equation  (2.4),  we  have 


where 


ur  =  {-uz)u  +  vuz2  , 


ikn(n  -  1) 

(-uz)  =  a(kQ,  r,z)  3  — - 


T 


v  3  bfk0-r>z)  3  w;  ■ 

Substituting  r  and  z  in  the  exact  solution,  we  have 

-1 


u(z,r)  3  \1  +  e 


z _ r_ 

2v  "  4v 
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The  initial  starting  field  and  surface  and  bottom  boundary  conditions  are 
then  determined  as  follows: 

Initial  Starting  Field 

u(z,r)  =  (: 

Surface  Condition 

u(o,r)  =  ^1 

Bottom  Condition 

u(Zmax,r)  = 

Solving  for  (-uz)  gives 


2 

where  n  = 

c  is  the  reference  sound  speed 
o 

c.  is  the  speed  of  sound  at  depth  z. 
Solving  for  ci  gives 


Because  of  the  constraints  placed  on  the  solution  of  Burger's  equation, 
this  test  example  has  no  real  significance  other  than  to  test  the  accuracy 


2z-r 


I  +  e 


ikft  V1 

(-2z+r)T 


irk. 


-I 


+  e 


<  ik 

—£■  (-2Zmax+r) 
,l*e2  / 


»Vl 
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of  the  IFO  model.  One  other  important  feature  of  this  test  example  is  that 
it  requires  the  user  to  supply  subroutines  UFIELD,  USVP,  BCON,  and  SCON. 

A  comparison  of  the  IFD  and  exact  solutions  at  approximately  7  meters 
in  range  is  shown  in  table  1. 


Table  1.  Comparison  of  IFD  and  Exact  Solutions 


I 

Z(I) 

U(I) 

IFD  10 
Exact 

0.10 

(0.49999E+00 
(0. 50000E+00 

-0.43163E+00) 

-0.43165E+00) 

20 

0.20 

(0.49998E+00 

(0.50000E+00 

-0.41366E+00) 

-0.41370E+00) 

30 

0.30 

(0.49998E+00 

(0.50000E+00 

-0.39631E+00) 

-0.39635E+00) 

40 

0.40 

(0.49997E+00 

(0.50000E+00 

-0.37953E+00) 

-0. 37958E+00) 

50 

0.50 

(0.49998E+00 

(0.50000E+00 

-0.36328E+00) 

-0.36333E+00) 

60 

0.60 

(0.49998E+00 

(0.50000E+00 

-0.34753E+00) 

-0.34756E+00) 

70 

0.70 

(0.49998E+00 
(0. 50000E+00 

-0. 33222E+00) 
-0.33225E+00) 

80 

0.80 

(0.49999E+00 

(0.50000E+00 

-0.31733E+00) 

-0.31736E+00) 

90 

0.90 

(0.49999E+00 
(0. 50000E+00 

-0. 30285E+00) 
-0.30286E+00) 

100 

1.00 

(0. 50000E+00 
(0.50000E+00 

-0.28872E+00) 

-0.28872E+00) 

The  input  ru.'stream  and  user-written  subroutines  UFIELD,  USVP,  BCON, 
and  SCON  are  listed  below. 

Input  Runstream  for  Exact  Solution 

100  .5  1500  1  0  1  100  0  1  1 
7.1  .001  0  .1  1  .1  0  0  0 
0  1 
100  1 
-1,-1 
0 
1 
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SUBROUTINE  UFIELO 


»  USER  STARTING  FIELD 

•  USER  WRITES  THIS  SUBROUTINE  IF  CAUSSIAN  FIELD  MOT  DESIRED 
»  UFIELD  IS  CALLED  IF  INPUT  PARAMETER  ISF  IS  HOT  ZERO 


•  UFIELD  SUBROUTINE  SUPPLIES: 

U  -  COMPLEX  STARTING  FIELD 


PARAMETER  MXLYR=101 ,  MXN  = 10000, MX SVPs 101 , MXTRKs 101 , NIUsI , 

C  MOUs  2  ,  HPUs 6 

COMPLEX  ACOFX, ACOFY.DCOF.BOTX, BOTY.BT A. HNK.HNKL.SU RX.SURY. TEMP, 

C  U.X.Y 

COMMON  /IFDCOM/ACOFX. ACOFY, ALPHA , BCOF . BETA( KXLYR) .BOTX.BOTY, 

C  BTA(HXN)  ,  CO.CSVPC  MXSVP)  ,  DR  .  DR  1  .  DZ  .  FR  <5 , 1 HN  K  .  ISF  .  I TY  PEB  . 

C  I  TYPES. IXSVP(MXLYR)  . KSVP . N . N 1 . NLYR . NSVP . NWSVP . R 1 2< MXN )  ,  R A  , 

C  RHO(MXLYR) , RSV P , SUR X . SUR Y , THETA , TRACK ( MXTR K , 2 )  .U(MXN) . 

C  X <  UXM ) .XKO.Y(MXN) , ZA . ZLYR ( MXLYR) , ZP, ZS , ZSVP( MXSVP) 

DATA  PI/j.^lS^fiSA/  ,  DEG/5  7 .2957'’./ 

C 

C  •••  STARTING  FIELD  FOR  EXACT  SOLUTION  TO  BURGER'S  PROBLEM 

DO  10  I  s  1  ,  N 
Z I = I *DZ 

U<  D  *  1 .0/(  1 . 0+CEXP(CMPLX{0 .0.  .5»XK0»(-2.0»ZI*RA) ) ) ) 

10  CONTINUE 
RETURN 
END 
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SUBROUTINE  USVP 

C  ****************************************************************** 

C  ***  USER  SOUND  VELOCITY  PROFILE  SUBROUTINE 

C  SUBROUTINE  USVP  IS  CALLED  EACH  DR  IN  RANGE  AS  LONG  AS 

C  KSVP  IS  NOT  ZERO.  KSVP  MAY  BE  USED  BY  USER  TO  TRANSFER  CONTROL 

C  IN  THIS  SUBROUTINE.  USER  INSERTS  LOGIC  TO  CLEAR  KSVP 

C  WHEN  USVP  IS  NO  LONGER  NEEDED.  IF  KSVP  NOT  CLEARED  BY  USER, 

C  USVP  IS  CALLED  EACH  STEP  IN  RANGE  UNTIL  RA  =  NEXT  RSVP. 

c  ****************************************************************** 

C  ***  USVP  SUBROUTINE  RETURNS: 

C  NLYR  -  NUMBER  OP  LAYERS.  LAYER  1  IS  WATER.  OTHERS  ARE  SEDIMENT 

C  ZLYR  -  ARRAY  -  DEPTH  OF  EACH  LAYER.  FIRST  IS  DEPTH  0?  WATER. 

C  RHO  -  ARRAY  -  DENSITY  OF  EACH  LAYER.  GRAMS/CUBIC  CM 

C  BETA  -  ARRAY  -  ATTENUATION  IN  EACH  LAYER.  DB/WAVELENGTH 

C  TXSVP  -  ARRAY  -  CONTAINS  POINTERS.  POINTS  TO  LAST  VALUE  OF  SVP 

C  IN  CORRESPONDING  LAYER.  SVP  TS  STORED  IN  ARRAYS  ZSVP 

C  AND  CSVP.  IXSVP(l)  POINTS  TO  LAST  SVP  POINT  IN  WATER. 

C  NSVP  -  NUMBER  OF  POINTS  IN  ZSVP  AND  CSVP.  ZSVP  AND  CSVP 

C  CONTAIN  THE  PROFILES  FOR  ALL  LAYERS. 

C  ZSVP  -  ARRAY  -  SVP  DEPTHS  -  METERS 

C  CSVP  -  ARRAY  -  SOUND  SPEED  -  METERS/SEC 

C  KSVP  -  AS  DESCRIBED  ABOVE. 

c  ****************************************************************** 

C 

PARAMETER  MXLYR5*  1 01,  MXN* 1 0000  ,  MXSVP*  101  ,MXTRK*101  ,NTU*1  , 

C  N0U=2,NPU*« 

COMPLEX  ACOFX , ACOFY , BCOF , BOTX , BOTY , BTA , HNK , HNKL , SURX , SURY , TEMP , 

C  U  ,  X ,  Y 

COMMON  /IFDCOM/ACOFX, ACOFY, ALPHA, BCOF, BETA(MXLYR) , BOTX, BOTY, 
BTA(MXN)  ,C0 .CSVP(MXSVP)  , DR , OR  1 , DZ , FRO , T HNK , ISF , TTYPEB , 
ITYPES, IXSVP( MXLYR) , KSVP , N , N 1 , NLYR , NSVP , NWSVP , R1 2 (MXN) , PA , 
RHO(MXLYR) , RSVP, SURX, SURY, THETA, TRACK (MXTRK , 2 ) ,U(MXN) , 

X ( MXN ) , XKO , Y ( MXN ) , ZA , ZLYR (MXLYR ) , ZP, ZS , ZSVP ( MXSVP) 

DATA  PT/3.1'1IS°2SS'1/, DEG /S7.ZQ5"’S/ 

C 

go  to  (inn, 2no, mn, 400)  ,ksvp 

nsvp*o 

RETURN 

inn  continue 
c 

C  ***  IF  KSVP*1,  CONTROL  TS  TRANSFERRED  HERE.  USER  LOADS 

C  NLYR.ZLYR(I) ,RHO(I) ,BETA(I) ,  AND  IXSVP(I)  WHERE  1*1 , NLYR. 

C  USER  ALSO  LOADS  NSVP , ZSVP ( I ) ,  AND  CSVP(T)  WHERE  1*1, NSVP. 

C  KSVP  may  BE  ALTERED  DEPENDING  ON  USER  LOGIC. 

C 

C  ***  SVP  FOR  EXACT  SOLUTION 

C 

NLYR= 1 
ZLYR(1)=1.0 
RHO ( 1 )  *  1 .0 
BETA ( 1 )  *0  1 
NSVP*101 

DZSVP=ZLYR(l)/(NSVP-l) 

XKP=?.n*PI*FRO/cn 
DO  lin  1*1, NSVP 


4-5 


TR  6659 


21 *  { X - 1 )  *DZSVP 

CSVP ( I ) *CO*SQRT (1 . 0+1 . 0/ ( COS (XKO* ( 2 1- . 5*RA) ) ) ) 
ZSVP ( I ) *2 1 
110  CONTINUE 

I XSVP { 1 ) *NSVP 
RETURN 
C 

200  CONTINUE 

C  ***  USER  INSERTS  CODE  HERE  IF  DESIRED 
RETURN 
C 

300  CONTINUE 

C  ***  USER  INSERTS  CODE  HERE  IF  DESIRED 

RETURN 
C 

400  CONTINUE 

C  ***  USER  INSERTS  CODE  HERE  IF  DESIRED 

RETURN 
END 


SUBROUTINE  8C0N 

C  •••  USER  PREPARED  BOTTOM  CONDITION  SUBROUTINE 

C  BCON  IS  CALLED  IF  INPUT  PARAMETER  ITYPEB  s  1 

C  SEE  MAIN  PROGRAM  FOR  DEFINITIONS 

C 

C  •••  SUBROUTINE  RETURNS: 

C  BOTY. BOTX 

c 

c 

PARAMETER  MXLYR* 1 0 1 . MXN* 1 0000 , MXSVPs 1 0 1 , MXTRK* 101 , NIUs  1 , 

C  NOU*  ? , NPU*6 

COMPLEX  ACOFX, ACOF Y . BCOF . BOTX , BOTY . BTA , HNK , HNKL, SURX , SURY, TEMP . 

C  U.X.Y 

COMMON  /IFDCOM/ ACOFX. ACOF Y. ALPHA, BCOF. BETA( KXLYR) .BOTX, BOTY, 

C  BTA( MXN ) , CO.CSVP( MXSVP) , DR , DR  1 , DZ.FRC, IHN K, ISF, ITYPEB, 

C  ITYPES, IXSVP(MXLYR) , KSVP , N , N 1 , NLYR , NSVP , NWS VP , R 1 2( MXN )  .  RA  , 

C  RHO(MXLYR) , RSVP, SURX, SURY, THETA, TRACK! MXTRK, 2) ,U(MXN) , 

C  X(MXtl)  .XKO.Y(KXN)  ,  ZA  ,  Z LYR(  MXLYR )  ,  ZP ,  ZS  .  ZSVP (  MXSVP) 

DATA  PI/ 3.  1U  159265*1/.  DEG/57. 29  57  A/ 

C 

IF(TKETA)  50.100,150 
C 

C  •••  THETA  LESS  THAN  0.0.  BOTTOM  SLOPES  UP. 

50  CONTINUE 

BOTYsU! N) 

C  BOTX  * . 

RETURN 

C  •••  THETA  :  0.0.  BOTTOM  IS  FLAT. 

C 

C  •••  BOTTOM  CONDITION  FOR  EXACT  SOLUTION  TO  BURGER’S  PROBLEM 

100  CONTINUE 

BOTYsU! N) 

BOTXsl  ,0/(  1  ,P*CF.XP(  CMPLX(  0.0,  .5"XK0»(-2.0»ZA*RA)))) 

RETURN 

C 

C  •••  THETA  GREATER  THAN  0.0,  BOTTOM  SLOPES  DOWN. 

150  CONTINUE 

C  BOTY*.... . 

C  SCTXs . 

RETURN 

END 
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SUBROUTINE  SCON 
C 

c 
c 
c 
c 
c 
c 

PARAMETER  MXLYRs 1 0 1 , MX  Ns  1 0000 , MXS VPs  1 0 1 , MXTRK* 1 0  1  , NI Us  1  . 

C  NOUs  2 , NPU  =  6 

COMPLEX  ACOFX.ACOFY.BCOF.BOTX.BOTY.BTA.HNK.HNKL, SURX. SURY. TEMP. 

C  U  .  X  .  Y 

COMMON  /IFDCOM/ACOFX, ACOFT . ALPHA , BCOF , BETA( MXLYR ) .BOTX.BOTY, 

C  BTA(MXN) , CO , CSVP( MXSVP) . DR , DR  1 , DZ , FR 0 . IHN K , ISF . ITY PEE , 

I  TYPES, IXSVP(MXLYR) ,KSVP, N. N1 , NL YR , NS V P . NWS VP , R 1 2( KXN) , RA. 
RHO(MXLYR) , RS VP . SUR X , S UR Y , THETA . TR ACK( MXTR K . 2) ,U(MXH)  , 
X(KXN) .XKO.Y(HXN) , 2A , ZLYRC HXLYR) . ZP . ZS , ZSVP( MXSVP) 

DATA  PI/ 3. 1M 159265 M/ , DEG/57.2957  8/ 

C 

IF( I  TYPES. NE.O)  GO  TO  100 
C 

C  •••  PRESSURE  RELEASE  SURFACE 

SUR  Y  =  0 . 0 
SURXsO.O 
RETURN 
C 

C  •••  USER  SURFACE  CONDITION 

C  •••  SURFACE  CONDITION  FOR  EXACT  SOLUTION  TO  BURGER'S  PROBLEM 

100  CONTINUE 

SURYsl .0/( 1 .0*CEXP( CHPLXIO.O,  .  5*XK0» ( R A-DR ) ) ) ) 

SU  RX  s  1 . 0/  (  1 ,0*CEXP(CMPLX(0.0.  ,5*XJC0#(  ♦RA)  ) )  ) 

RETURN 

END 


•••  SURFACE  CONDITION  SUBROUTINE 

IF  ITYPES  s  0,  SCON  SETS  SURY  AND  SURX  s  0.0. 

IF  ITYPES  NOT  0,  THE  USER  MUST  SUPPLY  SURY  AND  SURX. 
SEE  MAIN  PROGRAM  FOR  DEFINITIONS 
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4.2  RANGE-INDEPENDENT  PROBLEMS 

The  selection  of  range- independent  problems  includes  treatment  of  an 
isovelocity  shallow  water  environment  and  a  horizontal  interface. 


4.2.1  Isovelocity  Shallow  Water 

This  problem,  published  by  Jensen  and  Kuperman,14  considers  a  simple 
isovelocity  shallow  water  environment.  The  sound  speed  in  the  water  is 
1500  m/s.  The  water  depth  is  100  m,  and  both  source  and  receiver  are  placed 
50  m  in  depth.  In  the  bottom,  the  sound  speed  is  1550  m/s,  density  is 

1.2  g/cm3,  and  the  attenuation  is  1  dB/wavelength.  The  source  frequency  is 
500  Hz.  The  propagation  path  is  up  to  25  km  in  range. 

Solutions  obtained  by  Jensen  and  Kuperman,  using  a  normal  mode  model 
(SNAP)  and  a  parabolic  equation  split-step  model  (PAREQ)  developed  at 
SACLANT  Centre,  are  compared  with  the  solution  obtained  with  the  IFD  model. 
As  shown  in  figure  5,  all  solutions  are  in  excellent  agreement. 

The  input  IFD  runstream  that  produced  these  results  is  listed  below. 


Input  Runstream 

500  50  0  0  0  250  500  0  3  0 
25000  5  50  50  5000  50  0  0  0 
0  100 
25000  100 
-1,-1 
0 
0 
2 

100  1.0  -1.0 
0  1500 
100  1500 
200  1.2  1.0 
100  1550 

200  1550 

Note  that  although  the  bottom  parameters  were  extended  down  to  200  m,  the 
maximum  depth  of  the  solution  was  extended  to  250  m  as  requested  by  the 
combination  of  input  parameters  ZA  and  ITYPEB.  Artificial  attenuation 
was  then  applied  to  the  bottom-most  50  m  as  described  by  Brock.10 


4.2.2  Horizontal  Interface 


This  problem,  suggested  by  Dr.  H.  Bucker  of  NOSC  in  a  personal  com¬ 
munication,  considers  propagation  in  a  region  where  the  sound  speed  profile 
is  as  depicted  in  figure  6.  Source  and  receiver  depths  are  30  m  and  90  m, 
respectively.  Source  frequency  is  100  Hz.  At  240  m  in  depth,  the  density 
changes  abruptly  from  1.0  g/cm3  to  2.1  g/cm3.  No  attenuation  is  applied. 
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SOUND  SPEED  (m/s) 

1480  1490  1500  1510 


Figure  6.  Horizontal  Interface  Problem 

The  solution  produced  by  the  IFD  model  and  a  normal  mode  solution 
provided  by  Dr.  Bucker  are  compared  in  figure  7. 


Figure  7.  Solution  of  Horizortal  Interface  Problem 
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The  input  runstream  that  produced  these  results  is  listed  below.  It 
should  be  noted  that  the  bottom  was  artificially  extended  to  1200  m. 


Input  Runstream  for  Horizontal  Interface  Problem 

100  30  0  0  0  1200  600  030 
20000  2  50  90  10000  50  0  0  0 
0  240 
20000  240 
-1.-1 
0 
0 
2 

240  1.0  0.0 
0  1500 
120  1498 
240  1500 
512  2.1  0.0 
240  1505 
512  1505 


4.3  RANGE- DEPENDENT  PROBLEMS 

The  selection  of  range-dependent  problems  includes  treatment  of  depth 
dependent  environments,  interface  conditions,  and  a  homogeneous  Neumann 
bottom  boundary  condition. 


4.3.1  Shallow- to- Deep  Water  Problem 

This  problem,  extracted  from  Jensen  and  Kuperman,14  considers  propaga¬ 
tion  in  shallow-to-deep  water  as  shown  in  figure  8.  The  region  of  propaga¬ 
tion  is  bounded  by  a  pressure  release  surface  and  an  irregular  bottom  where 
the  bottom  remains  flat  at  50  m  in  depth  for  the  first  10  km;  at  10  km  the 
bottom  begins  to  slope  downward  until  it  levels  off  and  remains  flat  at  350  m 
in  depth.  Propagation  loss  was  calculated  at  two  different  sloping  angles, 
0.85  and  8.5  degrees.  The  sound  speed  in  the  water  is  1500  m/s;  in  the  bottom 
it  is  1600  m/s.  In  the  bottom,  a  density  of  1.5  g/cm3  and  an  attenuation  of 
0.2  dB/wavelength  are  used.  The  source  frequency  is  25  Hz,  and  the  source 
and  receiver  are  both  placed  at  a  depth  of  25  m. 

The  results  produced  by  SNAP,  PAREQ,  and  the  IFD  model  are  shown  in 
figures  9  and  10. 
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0  10  20  30  40  km 


Figure  8.  Shall ow-to-Deep  Water  Propagation 


The  input  runstream  that  produced  the  IFD  results  for  the  0.85  degree 
sloping  bottom  is  listed  below. 

Input  Runstream  (0.85  degree  slope) 


25  25  1500  0  0  1000  1000  030 
40000  10  100  25  10000  25  0  0  0 
0  50 

10000  50 
30000  350 
40000  350 
-1,-1 
0 
0 
2 

50  1.0  -1.0 
0  1500 
50  1500 
750  1.0  .2 
50  1600 
750  1600 
10000 
1 

30000 

0 

2 

350  1.0  -1.0 
0  1500 
350  1500 
750  1.0  .2 
350  1600 
750  1600 
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As  requested  in  the  input  runs tream,  user  subroutine  USVP  is  called  to 
supply  sound  speed  profiles  over  the  region  from  10  to  30  km  in  range.  Sub¬ 
routine  USVP  is  included  below. 


SUBROUTINE  USVP 

c 

C  •••  USER  SOUND  VELOCITY  PROFILE  SUBROUTINE 

C  SUBROUTINE  USVP  IS  CALLED  EACH  DR  IN  RANGE  AS  LONG  AS 

C  KSVP  IS  NOT  ZERO.  KSVP  HAY  BE  USED  BY  USER  TO  TRANSFER  CONTROL 

C  IN  THIS  SUBROUTINE.  USER  INSERTS  LOGIC  TO  CLEAR  KSVP 

C  WHEN  USVP  IS  NO  LONGER  NEEDED.  IF  KSVP  NOT  CLEARED  BY  USER. 

C  USVP  IS  CALLED  EACH  STEP  IN  RANGE  UNTIL  RA  *  NEXT  RSVP. 

C 

C  USVP  SUBROUTINE  RETURNS: 

C  NLYR  -  NUMBER  OF  LAYERS.  LAYER  1  IS  WATER.  OTHERS  ARE  SEDIMENT 

C  ZLYR  -  ARRAY  -  DEPTH  OF  EACH  LAYER.  FIRST  IS  DEPTH  OF  WATER. 

C  RHO  -  ARRAY  -  DENSITY  OF  EACH  LAYER.  GRAMS/CUBIC  CM 

C  BETA  -  ARRAY  -  ATTENUATION  IN  EACH  LAYER.  DB/W AVELENGTH 

C  IXSVP  -  ARRAY  -  CONTAINS  POINTERS.  POINTS  TO  LAST  VALUE  OF  SVP 

C  IN  CORRESPONDING  LAYER.  SVP  IS  STORED  IN  ARRAYS  ZSVP 

C  AND  CSVP.  IXSVP(I)  POINTS  TO  LAST  SVP  POINT  IN  WATER. 

C  NSVP  -  NUMBER  OF  POINTS  IN  ZSVP  AND  CSVP.  ZSVP  AND  CSVP 

C  CONTAIN  THE  PROFILES  FOR  ALL  LAYERS. 

C  ZSVP  -  ARRAY  -  SVP  DEPTHS  -  METERS 

C  CSVP  -  ARRAY  -  SOUND  SPEED  -  METERS/SEC 

C  KSVP  -  AS  DESCRIBED  ABOVE. 

c 

PARAMETER  MXLYR* 1 0 1 . MXN* 1 0000 . MXS VPs  1 0 1 . MXTRK* 1 0 1 . NIU* 1 . 

C  NOU*2,NPU=6 

COMPLEX  ACOFX, ACOFY.BCOF. BOTX , BOTY . BTA , HNK . HNKL , SURX . SUR Y . TEMP . 

C  U.X.Y 

COMMON  /IFDCOH/ACOFX. ACOFY, ALPHA, BCOF.BETAt MXLYR) .BOTX, BOTY, 

C  BTA(MXN) .CO.CSVP(MXSVP) , DR , DR  1 , DZ , FRO. IHNK , ISF , ITYPEB , 

C  ITYPES. IXSVP< MXLYR) , KSVP , N , N 1 , NLYR . NS VP . NWSVP . R 1 2( MXN )  ,  RA  , 

C  RHO( MXLYR) . RS V P, SUR X , SUR Y , THETA . TR ACKC MXTRK . 2) ,U(MXN) , 

C  X(MXN) .XKO.Y(MXN) , Z A , ZLYR( MXLYR) , ZP . ZS , ZSVP( MXS VP) 

DATA  P 1/3.  14159265  4/ , DEG/57.2957 8/ 

C 

GO  TO  ( 100.200, 300, «00)  .KSVP 
NSVPsO 
RETURN 
C 

100  CONTINUE 
C 

C  •••  IF  KSVP* 1 ,  CONTROL  IS  TRANSFERRED  HERE.  USER  LOADS 

C  NLYR.ZLYR(I) ,RHO(I) , BET A( I ) ,  AND  IXSVP(I)  WHERE  1*1, NLYR. 

C  USER  ALSO  LOADS  NSVP , ZSVP( I ) .  AND  CSVP(I)  WHERE  1*1, NSVP. 

C  KSVP  MAY  BE  ALTERED  DEPENDING  ON  USER  LOGIC. 

C 

C  •••  USER  SUPPLIES  SVP 

C 

NLY  R  *  2 

ZLYR( 1 ) *50.0*< R A- 10000.0) »T AN (THETA) 

RHO ( 1 )*  1 .0 
BETA(  1 ) *  -  1 . 0 
ZSVP( 1)*0.0 
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CSVP< 1 ) s 1 500 . 0 
ZSVP ( 2) s  ZLYR ( 1) 

CS VP( 2 ) a  1 500 . 0 
IXSVP(  1 ) *2 
ZLYR( 2 ) a750 . 0 
RHO  ( 2 ) s  1 . 0 
BETA( 2)  a  .  2 
2SVP(3)sZSVP(2) 

CSVP( 3 ) *  1 600 . 0 
ZS VP ( 4 ) a75 0 . 0 
CSVP( 4 ) *  1 600 . 0 
NSVPa4 
IXSVP(2)a4 
RETURN 

C 

200  CONTINUE 

C  USER  INSERTS  CODE  HERE  IF  DESIRED 

RETURN 

C 

300  CONTINUE 

C  •••  USER  INSERTS  CODE  HERE  IF  DESIRED 

RETURN 

C 

400  CONTINUE 

C  •••  USER  INSERTS  CODE  HERE  IF  DESIRED 

RETURN 
END 
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4.3.2  Deep-to-Shallow  Water  Problem 

This  problem,  also  extracted  from  Jensen  and  Kuperman,14  considers 
propagation  in  deep- to-sha How  water  as  shown  in  figure  11.  The  environment 
is  identical  to  that  of  the  previous  problem  except  that  the  shallow  and 
deep  portions  are  reversed  such  that  the  bottom  slopes  upward  in  the  direc¬ 
tion  of  propagation.  The  results  produced  by  SNAP,  PAREQ,  and  IFD  are 
shown  in  figures  12  and  13. 


0  10  20  30  40  km 


Figure  11.  Deep-to-Shallow  Water  Propagation 
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The  input  runstream  and  user  subroutine  USVP  which  produced  the  IFD 
results  for  the  8.5  degree  slope  are  listed  below. 

Input  Runstream  (8.5  degree  slope) 

25  25  0  0  0  1000  1000  0  3  0 
40000  10  100  25  10000  25  0  0  0 
0  350 
10000  350 
12000  50 
40000  50 
-1,-1 
0 
0 
2 

350  1.0  -1.0 
0  1500 
350  1500 
750  1.0  .2 
350  1600 
750  1600 
10000 
1 

12000 

0 

2 

50  1.0  -1.0 
0  1500 
50  1500 
750  1.0  .2 
50  1600 
750  1600 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 

100 

c 

c 

c 

c 


SUBROUTINE  USVP 


USER  SOUND  VELOCITY  PROFILE  SUBROUTINE 

SUBROUTINE  USVP  IS  CALLED  EACH  DR  IN  RANGE  AS  LONG  AS 

KSVP  IS  NOT  ZERO.  KSVP  MAY  DE  USED  BY  USER  TO  TRANSFER  CONTROL 

IN  THIS  SUBROUTINE.  USER  INSERTS  LOGIC  TO  CLEAR  KSVP 

WHEN  USVP  IS  NO  LONGER  NEEDED.  IF  KSVP  NOT  CLEARED  BY  USER. 

USVP  IS  CALLED  EACH  STEP  IN  RANGE  UNTIL  RA  s  NEXT  RSVP. 


•••  USVP  SUBROUTINE  RETURNS: 

NLYR  -  NUMBER  OF  LAYERS.  LAYER  1  IS  WATER.  OTHERS  ARE  SEDIMENT 
ZLYR  -  ARRAY  -  DEPTH  OF  EACH  LAYER.  FIRST  IS  DEPTH  OF  WATER. 
RHO  -  ARRAY  -  DENSITY  OF  EACH  LAYER.  GRAMS/CUBIC  CM 
BETA  -  ARRAY  -  ATTENUATION  IN  EACH  LAYER.  DB/WAVELENGTH 
IXSVP  -  ARRAY  -  CONTAINS  POINTERS.  POINTS  TO  LAST  VALUE  OF  SVP 
IN  CORRESPONDING  LAYER.  SVP  IS  STORED  IN  ARRAYS  ZSVP 
AND  CSVP.  IXSVP(I)  POINTS  TO  LAST  SVP  POINT  IN  WATER. 
NSVP  -  NUMBER  OF  POINTS  IN  ZSVP  AND  CSVP.  ZSVP  AND  CSVP 
CONTAIN  THE  PROFILES  FOR  ALL  LAYERS. 

ZSVP  -  ARRAY  -  SVP  DEPTHS  -  METERS 
CSVP  -  ARRAY  -  SOUND  SPEED  -  METERS/SEC 
KSVP  -  AS  DESCRIBED  ABOVE. 


PARAMETER  MXLYRs 1 0 1 . MXNs 1 0000 . HXS VPs  1 0 1 . MXTRKs 1 0 1 , N I  Us  1 . 

C  NOUs2,NPUs6 

COMPLEX  ACOFX. ACOF Y , BCOF , BOTX , BOTY . BTA . HNK , HNKL , SURX . SUR Y . TEMP . 

C  U.X.Y 

COMMON  /IFDCOM/ ACOFX. ACOFY, ALPHA. BCOF. BETA(MXLYR) .BOTX. BOTY. 

C  BTA(MXN) .CO.CSVP(MXSVP) . DR , DR  1 , DZ , FRO, IHNK. ISF . ITYPEB , 

C  ITYPES, IXSVP(MXLYR) . KSVP . N . N 1 , NLYR . NSVP , NVSVP . R 1 2( KXN ) , R A , 

C  RHO(MXLYR) ,RSVP,SURX.SURY.THETA,TRACK(MXTRK,2) ,U(MXN) . 

C  X(HXN) .XKO.Y(MXN) , ZA . ZLYRC MXLYR) , ZP . ZS , ZSVPC HXSVP) 

DATA  PI/  3.  I'M  5  9265*1/  .DEG/57.29578/ 

GO  TO  (  100.200.300.1100)  .KSVP 

NSVPsO 

RETURN 


CONTINUE 

•••  IF  KSVPsI,  CONTROL  IS  TRANSFERRED  HERE.  USER  LOADS 

NLYR.ZLYR(I) ,RHO(I) .BETA(I) .  AND  IXSVP(I)  WHERE  Isl.NLYR. 
USER  ALSO  LOADS  NS VP. ZSVP( I )  ,  AND  CSVP(I)  WHERE  Isl.NSVP. 
KSVP  MAY  BE  ALTERED  DEPENDING  ON  USER  LOGIC.  * 


•••  USER  SUPPLIES  SVP 


N  L  Y  R  s  2 

ZLYRC 1 )s350.0+( R A- 1 0000. 0)*TAN( THETA) 
RHO<  1  )s  1 .0 
BETA(  1 ) s-  1 . 0 
ZSVP(  DsO.O 
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CSVP< 1 )*1500.0 
ZSVP( 2) sZLYR (  1 ) 

CSVP( 2 ) s 1 500 . 0 
IXSVP( 1 ) =2 
ZLYR(2)*750.0 
RHO( 2) s  1 . 0 
BETA( 2) *  .  2 
ZSVP( 3) sZSVP( 2) 

CSVP{ 3 ) = 1 600 . 0 
ZSVP( 4)5750.0 
CSVP(  4 ) *  1 600 . 0 
NSVPs  4 
I XS VP( 2 ) s  4 
RETURN 

C 

200  CONTINUE 

C  USER  INSERTS  CODE  HERE  IF  DESIRED 

RETURN 

C 

300  CONTINUE 

C  •••  USER  INSERTS  CODE  HERE  IF  DESIRED 

RETURN 

C 

400  CONTINUE 

C  •••  USER  INSERTS  CODE  HERE  IF  DESIRED 

RETURN 
END 
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4.3.3  Shallow- to- Deep  Mater  Propagation  in  a  Wedge-Shaped  Region  With 
A  Rigid  Sloping  Bottom 

In  this  example,  the  IFD  model  was  used  to  propagate  the  acoustic  field 
in  a  wedge-shaped  region  with  a  rigid  sloping  bottom  as  shown  in  figure  14. 
The  purpose  of  this  example  is  to  exercise  the  optional  rigid  bottom  boundary 
condition  programmed  in  the  model.  For  this  case,  the  source  frequency  is 
80  Hz;  source  depth  is  15.2  m;  bottom  depth  at  the  location  of  the  source 
is  30.5  m;  and  the  sound  speed  profile  is  constant  at  1524  m/s.  The  bottom 
is  rigid  and  slopes  downward  at  an  angle  of  5  degrees.  The  initial  field 
propagated  by  the  IFD  model  was  generated  by  the  method  of  images15  and  is 
at  an  initial  range  of  348.4  m  from  the  source. 


STARTING  FIELD 
1 


Figure  14.  Shal low-to-Deep  Water  Propagation,  Wedge-Shaped 
Region  With  a  Rigid  Sloping  Bottom 


Numerical  results  were  compared  with  the  exact  solution  obtained  by  the 
method  of  images.  A  plot  of  propagation  loss  versus  range  at  a  receiver 
depth  of  27.4  m  is  given  in  figure  15. 
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Figure  15.  Propagation  Loss  Versus  Range,  Wedge-Shaped 
Region  With  a  Rigid  Sloping  Bottom 


The  initial  field  generated  by  the  method  images  consisted  of  400  points 
spaced  at  approximately  0.15  m  in  depth.  As  the  solution  was  marched  out 
in  range,  the  field  was  extended  deeper  and  deeper  until,  at  10  km,  the  field 
consisted  of  5740  points  in  depth. 

The  input  runstream  and  user  subroutine  UFIELD  which  produced  these 
results  are  included  below. 

Input  Runstream 

80  15.24003  0  1  348.3886  60.96012  400  1  0  0 
10000  10  50  .5  5000  27  0  0  0 
0  30.48006 
10000  905.38 
-1.-1 
0 
0 
1 

30.48006  1,0  0.0 
0  1524.003 
30.48006  1524.003 
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SUBROUTINE  UFIELD 


•••  USER  STARTING  FIELD 

•••  USER  WRITES  THIS  SUBROUTINE  IF  GAUSSIAN  FIELD  NOT  DESIRED 
•••  UFIELD  IS  CALLED  IF  INPUT  PARAMETER  ISF  IS  NOT  ZERO 


•••  UFIELD  SUBROUTINE  SUPPLIES: 

U  -  COMPLEX  STARTING  FIELD 


PARAMETER  NUU*3 

PARAMETER  MXLYR* 101 , MXN* 1 0000 . MXS VP* 1 0 1 , MXTRK* 1 0 1 . N IU« 1  , 

C  NOU*  2 , NPUs  6 

COMPLEX  ACOFX, ACOFY , BCOF , BOTX . BOTY . BTA , HNK , HNKL , SURX , SURY . TEMP , 

C  U.X.Y 

COMMON  /IFDCOM/ACOFX, ACOFY. ALPHA. BCOF,BETA( MXLYR) .BOTX.BOTY, 

C  BTA(MXN) , CO , CSVP( MXSVP) . DR , DR  1 , DZ . FRQ . IHN K , ISF . ITYPEB  , 

C  ITYPES.IXS VP (MXLYR) . KSVP , N . N 1 , NLYR . NSVP . NWSVP . R 1 2( MXN ) ,RA, 

C  RHO(  MXLYR ) , RSV P . SUR X . SURY . THETA . TRACKC MXTRK . 2)  ,U(MXN)  , 

C  X(MXN) . XKO, Y ( MXN ) .ZA.ZLYRf MXLYR) , ZP . ZS . ZSVP( MXSVP) 

DATA  PI/ 3. IN  1592654/ .DEG/57 .29578/ 

•••  STARTING  FIELD  GENERATED  BY  WEDGE  PROGRAM.  CONSTANT  SVP. 

MUST  BE  DIVIDED  BY  HANKEL  FUNCTION. 

SET  IHNK  s  1  IN  IFD  INPUT  RUNSTREAM. 

CALL  ASSIGN(NUU. 'WEDGE2.FLD' ) 

•••  BYPASS  WEDGE  DATA 

READ(NUU)  NANG.F.C1 .ZS" .ZSBB, RMIN.RMAX.DRR. ZM IN. ZMAX.DZZ, PHI 

***  READ  WEDGE  STARTING  FIELD 

READ(NUU)  NZ,R,(U(I) ,I*1,8X) 

•••  NZ  IS  NUMBER  OF  DEPTHS.  -  SHOULD  BE  EQUAL  TO  N 

•••  R  IS  RANGE  IN  FT.  -  SHOULD  BE  EQUAL  TO  RA  IN  METERS 

WEDGE  REFERENCES  TO  1  METER  BY  ADDING  -0.0 . 0»  ALOG 1 0(  3 . 280P  3 3 ) 

PROGRAM  WHICH  PLOTS  IFD  SOLUTION  SHOULD  DO  SAME. 

CALL  CLOSE(NUU) 

RETURN 

END 
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5.  CONCLUSIONS 

The  IFD  model  exhibits  significant  advantages  for  solving  the  parabolic 
wave  equation  because  of  its  generality,  useful  capabilities,  unconditional 
stability  of  the  method,  and  efficient  computation  of  the  wave  field.  In 
the  present  design,  the  program  package  is  basic  in  structure  and  is  flexible 
enough  for  easy  generalization  and  modification. 

The  capabilities  for  handling  surface  and  bottom  boundary  conditions  as 
well  as  horizontal  interface  conditions  are  important  features  which  enhance 
the  parabolic  equation  model.  Further  capabilities  can  be  incorporated  with¬ 
out  much  difficulty. 

The  solution  of  the  system  (equation  (2.21))  by  a  sparse  matrix  technique 
is  economical  and  efficient.  The  disadvantage  of  the  IFD  method  is  that  the 
boundary  condition  information  at  the  next  range-level  must  be  known.  In 
cases  where  the  bottom  condition  is  unknown,  the  option  to  extend  the  bottom 
with  an  artificial  absorbing  layer  such  that  the  field  at  the  bottom  becomes 
zero  may  be  helpful. 

This  report  presented  the  status  of  the  IFD  model  as  of  this  writing. 
Future  capabilities  to  be  built  into  the  model  include  wide  angle  propagation, 
multiple  irregular  interfaces,  automatic  step-size  determination,  high  fre¬ 
quency  propagation,  shear  waves,  and  others. 
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Appendix  A 
IFO  COMPUTER  LISTING 


A-l/A-2 
Reverse  Blank 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


•  IHPLICIT  FINITE-DIFFERENCE  METHOD  FOR  SOLVING  THE 

•  PARABOLIC  EQUATION  :  UR*A»U*B»UZZ  ;  WHERE 

•  A  =  I».5»KO*(N»tJ-1)  ;  AND  B=I*.5/KO 


•  D.  LEE  ANDG.  BOTSEAS,  CODE  3342 

•  NAVAL  UNDERWATER  SYSTEMS  CENTER 

•  NEW  LONDON,  CONNECTICUT  06320,  L'.S.A. 


VAX-11/780  VERSION  -  FORTRAN  IV* 


c 

*  ftft  « 

c 

*•• 

ACOFX 

c 

•  •• 

ACOFY 

c 

•  ft  • 

ALPHA 

c 

•  •• 

ATT 

c 

•  •• 

BCOF 

c 

•  •• 

BETA 

c 

**  • 

POTX 

c 

•  •• 

BOTY 

c 

•  •• 

BTA 

c 

•  •• 

CO 

c 

*•  * 

CSVP 

c 

ft  •• 

DEG 

c 

•  •ft 

DR 

c 

•  •• 

DR  1 

c 

•  ft  ft 

DZ 

c 

•  ft* 

DZZ 

c 

•  ft  ft 

FRQ 

c 

«  •• 

HNK 

c 

•  ft  ft 

HNKL 

c 

ft  ftft 

I BOT 

c 

c 

c 

•  ft  ft 

IDIAG 

c 

•  •• 

IHNK 

COEFFICIENT  'A'  AT  BOTTOM  -  AT  ADVANCED  RANGE 
COEFFICIENT  'A'  AT  BOTTOM  -  AT  PRESENT  RANGE 
VOLUME  ATTENUATION  -  DB/METER 

ATTENUATION  COEFFICIENT  FOR  ARTIFICIAL  ABSORBING  LAYER 

COEFFICIENT  'B*  -  RANGE  INDEPENDENT 

ARRAY  -  ATTENUATION  IN  LAYERS  -  DB/ WAVELENGTH 

COMPLEX  PRESSURE  AT  BOTTOM  AT  ADVANCED  RANGE  RA+DR 

COMPLEX  PRESSURE  AT  BOTTOM  AT  PRESENT  RANGE  RA 

ARRAY  -  PARTIAL  SOLUTION  OF  SYSTEM  OF  EQUATIONS 

REFERENCE  SPEED  OF  SOUND  -  METERS/SEC 

ARRAY  -  SOUND  VELOCITY  -  METERS/SEC 

CONVERSION  FACTOR  -  DEGREES/ RADIAN 

RANGE  STEP  -  METERS 

LAST  DR  USED  IN  ROUTINE  DIAG 

DEPTH  INCREMENT  OF  SOLUTION  -  METERS 

DEPTH  INCREMENT  FOR  ADJUSTING  LAYER  DEPTHS  IN  SLOPING  BOTTOM 
FREQUENCY  -  HZ 
HANKEL  FUNCTION  H0< 1 5 

EXTERNAL  FUNCTION  -  COMPUTES  HANKEL  FUNCTION  H0(1) 

BOTTOM  DEPTH  PRINT  FLAG 
s  0  -  DO  NOT  PRINT  BOTTOM  DEPTHS 
s  1  -  PRINT  BOTTOM  DEPTHS 
DIAGONAL  UPDATE  FLAG 
HANKEL  FUNCTION  FLAG 
=  0 


1 0* LOG ( R )  ADDED  TO 


•••  ILYR 

•••  IPZ 

ISF 


ISFLD 


HANKEL  FUNCTION  NOT  USED. 

SOLUTION . 

=  1  -  STARTING  FIELD  DIVIDED  BY  HANKEL  FUNCTION. 

SOLUTION  MULTIPLIED  BY  HANKEL  FUNCTION  BEFORE 
COMPUTING  PROPAGATION  LOSS. 

INDEX  FOR  ARRAYS  BETA,  ZLYR,  AND  RHO 
EVERY  IPZ'TH  VALUE  IN  DEPTH  IS  PRINTED 
STARTING  FIELD  FLAG 

*  0  -  PRO'.RAM  GENERATES  GAUSSIAN  STARTING  FIELD 
AT  RANGE  *  0.0.  SEE  SUBROUTINE  SFIELD. 

:  1  -  USER  SUPPLIES  STARTING  FIELD.  SEE  SUBROUTINE 
UFIELD. 

STARTING  FIELD  PRINT  FLAG 
:  C  -  DO  NOT  PRINT  STARTING  FIELD 
r  1  -  PRINT  STARTING  FIELD 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


•••  ISVP  - 


•••  ITEMP  - 
ITRK  - 
•••  ITYPEB- 
0  - 


2  -  NOT  RIGID 
2  -  NOT  RIGID 


•••  ITYPES- 


•*«  IWZ 

•••  IXSVP 


KSVP 

0 

1 

2 


SVP  PRINT  FLAG 

s  0  -  DO  NOT  PRINT  SOUND  VELOCITY  PROFILE 
:  1  -  PRINT  SOUND  VELOCITY  PROFILE 
TEMPORARY  VARIABLE 
INDEX  FOR  ARRAY  TRACK 
TYPE  OF  BOTTOM 

RIGID  BOTTOM  -  PROGRAM  SUPPLIES  BOTTOM  CONDITION 
-  NOT  RIGID  -  USER  SUPPLIES  BOTTOM  CONDITION 

-  SEE  SUBROUTINE  BCON 
NOT  RIGID  -  ABSORBING  LAYER  INTRODUCED 

-  FOLLOWS  CONTOUR  OF  BOTTOM 
NOT  RIGID  -  ABSORBING  LAYER  INTRODUCED  -  BUT 

-  BOTTOM  OF  ABSORBING  LAYER  KEPT  FLAT 
TYPE  OF  SURFACE 

:  0  -  PRESSURE  RELEASE.  SCON  SETS  SURY  AND  SURX  *  0.0 
*  .NE.  0  -  USER  INSERTS  CODE  IN  SCON  TO  COMPUTE  SURY  AND  S 
INDEX  INCREMENT  OF  RECEIVER  SOLUTIONS  TO  BE  WRITTEN  ON  DIS. 
ARRAY  OF  POINTERS  WHICH  POINT  TO  ENTRIES  IN  CSVP  AND  ZSVP 
IXSVP(I)  POINTS  TO  BOTTOM  DEPTH  AND  SPEED  IN  LAYER  1 
I XSVP ( 2 )  POINTS  TO  BOTTOM  DEPTH  AND  SPEED  IN  LAYER  2 
ETC. 

SVP  PROFILE  ''LAG 

PROFILE  IS  ON  CARDS  -  SEE  SUBROUTINE  SVP 
USER  SUPPLIES  PROFILE  1  -  SEE  SUBROUTINE  USVP 
USER  SUPPLIES  PROFILE  2  -  SEE  SUBROUTINE  USVP 


N  -  USER  SUPPLIES  PROF ILE  H  -  SEE  SUBROUTINE  USVP 

•••  MLYR  -  TEMPORARY  -  NUMBER  OF  LAYERS  INVOLVED  IN  SPECIFIC  CALCULATION 

MM  -  INDEX  -  M M* 1  POINTS  TO  FIRST  VALUE  OF  ARTIFICIAL  ABSORBING 
LAYER  IN  ARRAY  U 

•••  MXLYR  -  PARAMETER  -  MAXIMUM  NUMBER  OF  LAYERS 

-  MAX  DIMENSION  OF  ARRAYS  BET A . RHO . ZL YR . I XSV P 
MXH  -  PARAMETER  -  MAXIMUM  DIMENSION  OF  C,  X,  Y.  R12  AND  U  ARRAYS 

•••  MXSVP  -  PARAMETER  -  MAXIMUM  DIMENSION  OF  ARRAYS  CSVP  AND  ZSVP 

•••  MXTRK  -  PARAMETER  -  MAXIMUM  DIMENSION  OF  ARRAY  TRACK 

•••  N  -  NUMBER  OF  EOUI-SPACED  POINTS  III  U 

-  INCLUDES  BOTTOM  POINT  -  DOES  NOT  INCLUDE  SURFACE  POINT 
•••  N 1  -  DIAGONAL  ELEMENTS  NT  THRU  N  WILL  BE  COMPUTED 

•••  NA  -  NUMBER  OF  POINTS  IN  ABSORBING  LAYER 

•••  NIU  -  PARAMETER  -  INPUT  UNIT  NUMBER 

MLYR  -  NUMBER  OF  LAYERS 

NOLD  -  NUMBER  OF  RECEIVER  DEPTHS  ON  ENTRY  TO  ROUTINE  CRNK 
•**  NOU  -  PARAMETER  -  OUTPUT  UNIT  NUMBER 
NPU  -  PARAMETER  -  PRINTER  UNIT  NUMBER 

•**  NSVP  -  NUMBER  OF  POINTS  IN  CSVP  AND  ZSVP 

•»*  HWSVP  -  NUMBER  OF  POINTS  IN  LAYER  t  SVP 

»»»  NZ  -  NUMBER  OF  SOLUTION  DEPTHS  TO  BE  WRITTEN  ON  DISK 

»•«  OLDR  -  RANGE  INCREMENT  AT  START  OF  PROBLEM 

•••  PDR  -  RANGE  INCREMENT  AT  WHICH  SOLUTION  IS  PRINTED  -  METERS 

•••  PDZ  -  DEPTH  INCREMENT  AT  WHICH  SOLUTION  IS  PRINTED  -  METERS 

PI  -  THE  VALUE  OF  PI 

•••  PL  -  PROPAGATION  LOSS  -  DB 

•»»  R 1  -  RANGE  AT  WHICH  BOTTOM  DEPTH  IS  AVAILABLE  -  METERS 

R 2  -  NEXT  RANGE  AT  WHICH  BOTTOM  DEPTH  IS  AVAILABLE  -  METERS 
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i 


C  •••  R  1  2 

C  •••  RA 

C 

C 

C  •••  R  A  +  DR 

C  •••  RHO 

C  •••  RMAX 

C  •••  RSVP 

C  •••  SURX 

C  SURY 

C  TEMP 

C  •••  THETA 

C 
C 
C 

C  •••  TM 

C  •••  TRACK 

C  •••  U 

C  •••  WDR 

C  •••  HDZ 

C 

c 

C  •••  X 

C  •••  XKO 

C  •••  XPR 

C  •••  XWR 

C  •••  Y 

C  Z1 

C  Z2 

C  •••  ZA 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

C  •••  ZLYR 

C  •••  ZI 

C  •••  ZS 

C  •••  ZSVP 


ARRAY  OF  DENSITY  RATIOS 

HORIZONTAL  RANGE  OF  STARTING  FIELD  FROM  SOURCE  -  METERS 
RA  IS  SET  TO  0.0  IF  STARTING  FIELD  IS  GAUSSIAN.  RA  IS 
INCREMENTED  AS  SOLUTION  IS  MARCHED  OUT  IN  RANGE. 

RANGE  TO  WHICH  SOLUTION  IS  TO  BE  ADVANCED  -  METERS 

ARRAY  -  DENSITY  IN  LAYER 

MAXIMUM  RANGE  OF  SOLUTION  -  METERS 

NEXT  RANGE  AT  WHICH  SVP  IS  AVAILABLE  -  METERS 

COMPLEX  PRESSURE  AT  SURFACE  AT  ADVANCED  RANGE  RA+DR 

COMPLEX  PRESSURE  AT  SURFACE  AT  PRESENT  RANGE  RA 

TEMPORARY  VARIABLE  -  COMPLEX 

SLOPE  OF  BOTTOM  -  RADIANS 

. E Q . 0  --  FLAT  BOTTOM 

. GT. 0  —  SHALLOW  TO  DEEP 

. LT . 0  —  DEEP  TO  SHALLOW 

ARRAY  -  TIME  OF  DAY 

2  DIM.  ARRAY  -  RANGE  AND  DEPTH  OF  WATER  -  METERS 
ARRAY  -  COMPLEX  ACOUSTIC  PRESSURE  FIELD 

RANGE  STEP  AT  WHICH  SOLUTION  IS  WRITTEN  ON  DISK  -  METERS 
DEPTH  INCREMENT  AT  WHICH  SOLUTION  IS  WRITTEN  ON  DISK  -  METERS 
WDZ  SHOULD  BE  SELECTED  SO  THAT  PLOT  PROGRAM  DOES  NOT 
INTERPOLATE  BETWEEN  WIDELY  SPACED  RECEIVERS. 

ARRAY  -  MAIN  DIAGONAL  OF  MATRIX  AT  ADVANCED  RANGE 
REFERENCE  WAVE  NUMBER 

RANGE  AT  WHICH  SOLUTION  IS  PRINTED  -  METERS 

RANGE  AT  WHICH  SOLUTION  IS  WRITTEN  ON  DISK  -  METERS 

ARRAY  -  MAIN  DIAGONAL  OF  MATRIX  AT  PRESENT  RANGE 

DEPTH  OF  WATER  AT  RANGE  R1  -  METERS 

DEPTH  OF  WATER  AT  RANGE  R?  -  METERS 

DEPTH  OF  FIELD  AT  RANGE  RA  -  METERS 

INITIAL  DEPTH  OF  STARTING  FIELD  AT  RANGE  RA  IS 

AS  FOLLOWS: 

IF  ITYPEB  *  0  OR  1.  ZA  IS  MAXIMUM  DEPTH  OF 
BOTTOM-MOST  SEDIMENT  LAYER  AT  INITIAL  RANGE  OF 
STARTING  FIELD.  IF  ITYPEB  :  2  OR  3,  ZA  IS  MAXIMUM 
DEPTH  OF  ARTIFICIAL  ABSORBING  LAYER  AT  INITIAL 
RANGE  OF  STARTING  FIELD.  PROGRAM  INSERTS  LAYER. 

RHO  AND  BETA  ARE  OBTAINED  FROM  LAYER  ABOVE. 

SPEED  IS  BOTTOM-MOST  SPEED  FROM  LAYER  ABOVE. 

AS  SOLUTION  PROGRESSES, 

ZA  IS  UPDATED  IF  OCEAN  BOTTOM  HOT  FLAT.  IF  ITYPEB  s  3. 

BOTTOM  OF  ABSORBING  LAYER  REMAINS  FLAT. 

ARRAY  -  DEPTH  TO  BOTTOM  OF  LAYER  -  METERS 
DEPTH  OF  RECEIVER  'I'  -  METERS 
SOURCE  DEPTH  -  METERS 

ARRAY  -  DEPTH  OF  SOUND  VELOCITY  -  METERS 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

INPUT 

c 

Mil 

c 

INPUT 

UNIT  NUMBER 

s  NIU 

c 

INPUT 

FILE  NAME 

=  IFD. 

IN 

c 

CONTENTS 

CARD  IMAGES  IN  FREE 

c 

CARD 

1 

FRQ.ZS, 

CO. ISF 

. RA , ZA, 

c 

CARD 

2 

RMAX, DR 

,  WDR , WDZ , PDR  , 

c 

CARD 

3 

R  1  ,  Z  1 

•  * 

c 

CARD 

u 

R2.Z2 

ft 

BOTTOM 

c 

# 

ft  ft 

RANGE  , 

c 

ft 

c 

CARD 

N 

, 

ft  ft 

c 

CARD 

M+1 

-1.-1 

c 

CARD 

N*2 

RSVP 

c 

CARD 

N+3 

KSVP 

c 

CARD 

N*4 

NLYR 

c 

CARD 

N*5 

ZLYR( I) 

. RHO(  I 

) . BETA( 

c 

CARD 

N*6 

ZSVP( 1 ) 

. CSVP( 

1  ) 

c 

CARD 

N*7 

ZSVP( 2) 

,CSVP( 

2) 

c 

. 

c 

• 

, 

c 

. 

c 

CARD 

N*M 

ZSVP( J) 

. CSVP( J ) 

ft 

ft  ft 

REPEAT 

ft  ft 

REPEAT 

ft*  ft 

FOR  EACH 

•  ft  ft 

FOR  EACH 

•  ft 

PROFILE 

ft  ft 

LAYER 

ft 

ft 

ft 

QUICK  REFERENCE  AND  MOTES  FOR  CARD  INPUT 
•»»  UNITS:  METERS  AND  METERS/SEC  EXCEPT  AS  NOTED 


FRO 

ZS 

CO 

ISF 

RA 

ZA 


IF  CO  =  O.O,  CO  IS  SET  TO  AVERAGE 


1  =  USER  FIELD. 


FREOUENCY  (HZ) 

SOURCE  DEPTH 
REFERENCE  SOUND  SPEED. 

SPEED  IN  FIRST  LAYER. 

STARTING  FIELD  FLAG.  0  =  GAUSSIAN. 

IF  ISF  s  0,  RA  IS  SET  TO  ZERO. 

HORIZONTAL  RANGE  FROM  SOURCE  TO  STARTING  FIELD. 

RA  IS  SET  TO  0.0  IF  ISF  =  0. 

DEPTH  OF  STARTING  FIELD  AT  RANGE  RA.  IF  ZA  :  0.0. 

TO  MAX  DEPTH  OF  BOTTOM  LAYER  IN  FIRST  PROFILE.  IF 
2  OR  3  AND  ZA  s  0.0,  ZA  IS  SET  TO  («/3)»MAX  DEPTH  OF  BOTTOM 
LAYER.  IF  ITYPEB  s  2  OR  3  AND  ZA  NOT  ZERO,  THE  ARTIFICIAL 
BOTTOM  LAYER  IS  EXTENDED  TO  ZA  METERS  PROVIDED  THAT  ZA 
IS  GREATER  THAN  OR  EQUAL  TO  MAX  DEPTH  OF  BOTTOM  LAYER 
III  FIRST  PROFILE. 

NUMBER  OF  EQUISPACED  RECEIVERS  IN  STARTING  FIELD.  IF  N  s  0, 
N  IS  SET  SO  THAT  THE  RECEIVER  DEPTH  INCREMENT  IS  LESS  THAN 


ZA  IS  SET 
ITYPEB  = 


IF  N  IS  GREATER  THAN  MXN , 


IHNK  r 


OR  EQUAL  TO  1/U  WAVELENGTH. 

N  IS  SET  TO  MXN. 

HANKEL  FUNCTION  FLAG.  IHNK  s  0,  DON'T  USE  HAMKEL  FUNCTION. 
IHNK  s  1 .  DIVIDE  STARTINC 
MULTIPLY  THE  SOLUTION  FII 


COMPUTING  PROPAGATION 
IHNK  SHOULD  3E  SFT  TO 


LOSS  . 
0.  IF 


0, 

DO 

N  '  T 

USE 

HAN 

ELD 

BY 

HA 

NKEL 

FUN 

BY 

HAN 

KEL 

FUNCTIO 

IF 

STA 

RTI 

NG  FT 

ELD 

STA 

RTI 

NG 

FIELD 

IS 
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I 


i 


i 


i 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


IHNK  SHOULD  BE  SET  TO  1. 

ITYPEB  s  TYPE  OF  BOTTOM  PROCESSING. 

s  0  -  RIGID  BOTTOM.  PROGRAM  SUPPLIES  BOTTOM  CONDITION. 

=  1  -  USER  SUPPLIES  BOTTOM  CONDITION.  USER  WRITES  SUBROUTINE 
BCON. 

=  2  -  ARTIFICIAL  ABSORBING  BOTTOM  INTRODUCED.  FOLLOWS  CONTOUR 
OF  WATER/BOTTOM  INTERFACE. 

s  3  -  ARTIFICIAL  ABSORBING  BOTTOM  INTRODUCED.  BOTTOM  OF 
LAYER  KEPT  FLAT. 

ITYPES  =  TYPE  OF  SURFACE 

s  0  -  PRESSURE  RELEASE.  SCON  SETS  SURY  AND  SURX  s  0.0 
NOT  0  -  USER  INSERTS  CODE  IN  SCON  TO  COMPUTE  SURY  AND  SURX 


RHAX  =  MAXIMUM  RANGE  OF  SOLUTION 

DR  =  RANGE  STEP.  IF  DR  =  0,  DR  IS  SET  TO  1/2  WAVELENGTH. 

IF  BOTTOH  OF  PROBLEM  IS  NOT  FLAT.  DR  IS  RECOMPUTED 
SO  THAT  MAX  DE/TH  IS  EITHER  INCREMENTED  OR  DECREMENTED 
BY  DZ.  SOLUTION  IS  COMPUTED  EVERY  DR  METERS. 

WDR  =  RANGE  STEP  AT  WHICH  SOLUTION  IS  WRITTEN  ON  DISK.  IF  WDR  NOT  0. 
AN  OUTPUT  DISK  FILE  IS  ASSIGNED. 

ROUNDED  TO  NEAREST  DR. 

W DZ  =  DEPTH  INCREMENT  AT  WHICH  SOLUTION  IS  WRITTEN  ON  DISK. 

ROUNDED  TO  NEAREST  DZ . 

PDR  =  RANGE  STEP  AT  WHICH  SOLUTION  IS  PRINTED. 

ROUNDED  TO  NEAREST  DR. 

PDZ  =  DEPTH  INCREMENT  AT  WHICH  SOLUTION  IS  PRINTED. 

ROUNDED  TO  NEAREST  DZ. 

ISFLD  s  0  -  DON'T  PRINT  STARTING  FIELD. 

=  1  -  PRINT  STARTING  FIELD. 

ISVP  s  0  -  DON'T  PRINT  SOUND  VELOCITY  PROFILE. 

=  1  -  PRINT  SOUND  VELOCITY  PROFILE. 

IBOT  =  0  -  DON'T  PRINT  BOTTOM  DEPTHS. 

=  1  -  PRINT  BOTTOM  DEPTHS. 


R1.Z1 
R  2 ,  Z2 
RSVP 
KS  VP 


NLYR 

ZLYR( I ) 
RHO ( I) 
BETA( I) 

ZSVP 
CSVP 
ZSVP(  1  ) 
CSVP(  1) 


=  BOTTOM  PROFILE.  FIRST  RANCE  AND  DEPTH  OF  WATER, 
s  ETC.  (-1,-1)  MARKS  THE  END  OF  THE  BOTTOM  PROFILE, 
s  RANGE  OF  SVP 
s  SVP  FLAG. 

s  0  -  SVP  IN  INPUT  RUNSTREAM 

s  NOT  ZERO  -  PROFILE  (CARDS  N*4  THRU  N-t-M)  IS  SUPPLIED  BY 
USER.  USER  WRITES  SUBROUTINE  USVP.  KSVP  MAY  BE  USED  IN 
COMPUTED  GOTO  STATEMENT  TO  TRANSFER  CONTROL  IN  USVP. 
s  NUMBER  OF  LAYERS.  IF  ITYPEB  :  2  OR  3.  PROGRAM  INSERTS 
AN  ARTIFICIAL  LAYER  AND  INCREMENTS  NLYR  BY  1. 
s  MAX  DEPTH  OF  LAYER  I  IN  PROFILE 
=  DENSITY  IN  LAYER  I  (G/CM»»3) 

=  ATTENUATION  IN  LAYER  I  ( DB/ WAVELENGTH ) 

IF  BETA(I)  IS  NEGATIVE,  ATTENUATION  IS  COMPUTED, 
s  DEPTH  ARRAY  FOR  SOUND  SPEED 
s  SPEED  ARRAY  FOR  SOUND  SPEED 
s  DEPTH  OF  WATER  TO  TOP  OF  LAYER  I 
s  SPEED  OF  SOUND  AT  TOP  OF  LAYER  I 


L  I 


J 
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C  ZSVP(J)  3  DEPTH  OF  WATER  TO  BOTTOM  OF  LAYER  I 

C  CSVP(J)  =  SPEED  OF  SOUND  AT  BOTTOM  OF  LAYER  I 

C  IF  ONLY  ONE  SVP  INPUTTED,  IT  IS  USED  THRU  ENTIRE  PROBLEM. 

C  IF  MORE  THAN  ONE  SVP  INPUTTED,  LAST  SVP  IS  USED  THRU  REMAINDER 

C  OF  PROBLEM. 

C 

c  ••••••■*•>•■••••••••••••••■••••••••••••••••••••••••••••••••••••••» 

C  •**  OUTPUT 

C  OUTPUT  UNIT  NUMBER  s  NOU 

C  OUTPUT  FILE  NAME  s  IFD.OUT 

C  CONTENTS:  AS  FOLLOWS  -  UNFORMATTED 

C 

C  FRO.ZS.CO. ISF,RA,ZA,N, IHNK, ITYPEB , ITYPES , RM AX , DR , WDR , DZ . NLYR , ZLYR . 

C  RHO , BETA 

C  NZ,RA,WDZ,(U(I) .IsIWZ.N.IWZ)  -  (FOR  EACH  WRITE  RANGE  WDR) 

C 

C  PRINTER  UNIT  NUMBER  s  NPU 

c 

PARAMETER  MXLYRs 1 0 1 . MXNs 1 OOOO . MXS VPs  1 0 1 , MXTRKs 1 0 1 . N IUs 1 , 

C  NOUs2,NPUs6 

COMPLEX  ACOFX. ACOFY. BCOF.BOTX, BOTY, BTA ,HNK, HNKL.SURX, SUR Y, TEMP, 

C  U.X.Y 

COMMON  /IFDCOM/ ACOFX, ACOFY. ALPHA, BCOF, BET A( MXLYR) .BOTX.BCTY, 

BTA(MXN) .CO.CSVP(MXSVP)  . DR , DR  1 . DZ , FR Q . IHN K , ISF , ITYPEB , 

ITYPES. IXS VP (MXLYR) . KS VP . N , N 1 . NLYR , NSVP . NWSVP , R 1 2( MXN ) , R A , 
RHO(MXLYR)  . RS V P . SUR X , SUR Y , THETA , TR ACK( MXTR K . 2)  ,U(MXH)  , 

X(MXN) .XKO.Y(MXN) . Z A , ZLYR ( MXLYR ) . ZP , ZS . ZSVP ( MXS VP ) 

DATA  PI/ 3. 1 41 592654/ , DEG/ 5 7. 29578/ 

BYTE  TM( 8 ) 

C 

C  •••  TIME  OF  DAY  AT  START  OF  RUN 

CALL  TIME(TM) 

10  FORMAT(5X.8A1/) 

C 

C  •**  READ  INPUT  PARAMETERS 

CALL  ASSIGN(NIU, ' IFD.IN' ) 

READ(NIU.»)  FRO.ZS.CO, ISF, R A, Z A, N, IHNK, ITYPEB,  'TYPES 
READ( NIU  ,»)  RHAX,DR,WDR,WDZ, PDR , PDZ , ISFLD , ISVP , IBOT 

C 

C  IF  GAUSSIAN  STARTING  FIELD,  RA  MUST  BE  0. 

IF( ISF.EQ.O)  RAsO.O 
C 

C  •••  READ  BOTTOM  PROFILE  -  RANGE, DEPTH 

DO  22  I s  1  , MXTR K 

8  £AD( N IU , * )  TRACK(I.I) .TRACK (1, 2) 

ITRKsI 

C  •••  END  OF  PROFILE? 

IF(TRACK( I, 1 ) . LT.O.O)  GO  TO  2’ 

C  •••  NO 

22  CONTINUE 
C  •••  ERROR? 
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23  IF( ITRK . £0. I.OR.ITRK.GT.MXTRK)  CO  TO  28 

C  NO.  EXTEND  LAST  DEPTH  BEYOND  MAX  RANGE. 

TRACK( ITRK. 1)s1 .OE*38 

TRACK ( ITRK.2) sTR ACK ( ITRK- 1 . 2) 

ITRK=  1 

R2=TR ACKC ITRK.  1 ) 

Z2  sTR  ACK ( ITRK.2) 

C  FIND  BOTTOM  SEGMENT  WHICH  CONTAINS  STARTING  RANGE 

25  R 1 s  R2 

Z  1  =Z2 

ITRKs ITRK+ 1 

I F ( ITRK.GT.MXTRK)  GO  TO  28 
R2  =  TR ACKC ITRK.  1 ) 

Z2sTRACK( ITRK.2) 

C  ADVANCE  TRACK  IF  NECESSARY  SO  THAT  R 1 . LE . R A . LT . R2 

IF(RA.GE.R2)  GO  TO  25 
C  R 1  MUST  BE  LE  RA 

IFCR1.LE.RA)  GO  TO  30 
WRITEC  N  PU , 27) 

27  FORMATC 1 X. ' DEPTH  OF  BOTTOM  AT  INITIAL  RANGE  MISSING') 

STOP 

28  CONTINUE 
ITEMPrMXTRK 
WRITECNPU.29)  ITEMP 

29  FORMATC IX, 'ERROR.  BOTTOM  PROFILE  MISSING  OR  TOO  MANY  POINTS. 

CM AX  IS  • .15) 

STOP 

30  CONTINUE 
C 

C  COMPUTE  SLOPE  OF  BOTTOM 

THETAsATAN2(Z2-Z1,R2-Rt) 

C 

C  •••  READ  RANGE  OF  SVP 

32  READ(NIU.».END=33)  RSVP 

C  •••  SVP  BEYOND  START  RANGE? 

IF(RSVP.GT.RA)  GO  TO  33 

C  NO.  DETERMINE  IF  SVP  IN  RUNSTREAM  OR  SUPPLIED  BY  SUBROUTINE  USVP. 

READ(NIU.*,ENDs33)  KSVP 

IF(KSVP.EG.O)CALL  S VP ( N LYR , ZL YR , RHO , BETA , IXS VP . NS VP , ZS VP , CS VP) 
IF(KSVP.NE.O)CALL  USVP 
C  •••  ERROR  IN  SVP? 

IF(NSVP.NE.O)  GO  TO  35 
C  •••  YES. 

33  WRITEC  NPU,  3<»)  RA 

34  FORMATC IX. ' SVP  MISSING  OR  INPUT  ERROR.  RANGE  s  '.F8.1,'  M.') 

STOP 

C 

35  CONTINUE 

C  •••  SAVE  POINTER  TO  LAST  SVP  VALUE  IN  FIRST  LAYER. 

NWSVPsIXSVPC 1) 

C 

C  •••  IF  CO  NOT  SPECIFIED.  SET  CP  TO  AVERAGE  SPEED  OF  FIRST  PROFILE 

C  IN  FIRST  LAYER  AT  INITIAL  RANGE 

IFCCO.NE.O.O)  GO  TO  4? 
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DO  MO  1:2 , MWSVP 

C0sC0*( ZSVP( I)-ZSVP( 1-1 ) )•( CSVP( 1-1 )*.5»( CSVP( I)-CSVP( 1-1 ) ) ) 

MO  CONTINUE 

C0=C0/ZSVP(NWSVP) 

M5  CONTINUE 

•••  COMPUTE  REFERENCE  WAVE  NUMBER 
XK0*2.0»PI»FR<3/C0 

COMPUTE  ATTENUATION  -  SACLANT  MEMO  SM-121  (JENSEN  +  FERLA) 

•••  MODIFIED  AS  FOLLOWS: 

•••  IF  INPUTTED  BETA  IS  LT  0.0,  ALPHA  IS  COMPUTED  IN  DB/METER 
AND  USED  FOR  BETA 

ALPHA=FRO»FRO»( .007*( . 155*1. 7)/( 1 .7*1 . 7*F RC»FRQ« . OOOOO 1 ) ) • 1 . 0E-09 

•••  ADJUST  LAYER  DEPTHS  IN  CASE  BOTTOM  SLOPES  AND  RA.NE.RSVP 
•••  ASSUMES  LAYERS  ARE  PARALLEL  AND  FOLLOW  BOTTOM  CONTOUR. 

•••  LAYER  DEPTHS  ARE  ENTERED  WITH  SVP. 

DO  50  ILYRs 1 , NL YR 

0  ZLYR( ILYR)sZLYR( I LYR ) ♦( RA-RSVP) »TAN( THETA) 

•••  GET  RANGE  OF  NEXT  SVP 
READ(NIU,».ENDs55)  RSVP 
GO  TO  56 

•••  ONLY  ONE  SVP  -  SET  RSVP  LARGE  SO  SAME  SVP  USED  FOR  WHOLE  PROBLEM 

5  RSVPs 1 . OE+38 

6  CONTINUE 

•••  IF  STARTING  FIELD  IS  BEYOND  NEXT  SVP.  GO  BACK  AND  GET  NEXT  SVP 
IF( RSVP.LE.RA)  GO  TO  32 

•••  IF  ITYPEB  =  2  OR  3,  AND  ZA  *  0,  EXTEND  BOTTOM  BY  SETTING  ZA  TO 
•••  M/3  MAX  DEPTH 

•••  SEE  NORDA  TECH  NOTE  12,  JAN  78.  H.  K.  BROCK 

ip  ITYPEB  s  2  OR  3  AND  ZA  NOT  0.  EXTEND  BOTTOM  TO  ZA  METERS 
IF( ITYPEB.LT. 2)  GO  TO  60 
I F ( ITYPEB. GT. 3)  GO  TO  60 

•••  EXTEND  BOTTOM  TO  ZA  METERS  IF  ZA  NOT  ZERO 
•••  EXTEND  BOTTOM  M/3  MAX  DEPTH  IF  ZA  s  ZERO 
IF( ZA.EQ.O.O)  ZA=M.0»ZLYR(NLYR)/3.0 
IF( ZA.GE.ZLYR( NLYR) )  GO  TO  53 

•••  USER  ATTEMPTED  TO  EXTEND  DEPTH  IN  NEGATIVE  DIRECTION 
WRITE(NPU,57) 

7  FORMAT? IX. ' ERROR.  ZA  RESET  TO  MAX  DEPTH  OF  BOTTOM  LAYER.') 

ZA=ZLYR( NLYR) 

•••  INSERT  PARAMETERS  FOR  EXTENDED  BOTTOM  IN  APPROPRIATE  ARRAYS 

8  NLYRsNLYR+1 

•••  STORE  DEPTH  OF  ARTIFICIAL  ABSORBING  LAYER 
ZLYR( NLYR) sZA 

C  •••  USE  SAME  DENSITY  AND  ATTENUATION  AS  LAYER  AEOVE 

RHC ( NLYR) sRHO( NLYR-1 ) 

BETA( NLYR) =BETA( NLYR-1 ) 

C  •••  USE  BOTTOM  SPEED  OF  LAYER  ABOVE 

NSVPsNSVP+1 
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IXSVP<  NLYR) =  NSVP 
ZSVP(NSVP)*ZLYR(NLYR) 

CSVP(NSVP)=CSVP(NSVP-1) 

GO  TO  62 
60  CONTINUE 

C  IF  BOTTOM  NOT  EXTENDED  AND  ZA  =  0,  SET  ZA  TO  MAX  DEPTH  INPUTTED 

C  WITH  FIRST  SVP 

IFC ZA.EQ.O.O)  ZAsZLYRC NLYR) 

C 

C  •••  IF  N  NOT  SPECIFIED.  SET  N  SO  THAT  DZ  .LE.  1/4  WAVELENGTH 

62  IFCN.EQ.O)  Ns(4»ZA«FRQ/C0)+1 
IF(N.LE.MXN)  GO  TO  65 
ITEMPsMXN 

WRITEC  NPU , 64)  ITEMP 

64  FORMATC '  ERROR.  N  TOO  LARGE.  N  RESET  TO  '.14*.') 

NsMXN 

65  CONTINUE 
C 

C  •••  COMPUTE  RECEIVER  DEPTH  INCREMENT  -  DZ  MAY  BE  SUCH  THAT  RECEIVERS  DO 

C  •••  NOT  LIE  EXACTLY  ON  LAYER  INTERFACES 

DZsZA/N 
C 

C  •••  IF  BOTTOM  IS  FLAT,  RANGE  STEP  MAY  BE  SPECIFIED 

C  •••  IF  RANGE  STEP  NOT  SPECIFIED.  SET  IT  EQUAL  TO  1/2  WAVELENGTH  !??? 

IF(DR.EQ.O.O)  DRs . 5*C0/FRQ 
C  IF  BOTTOM  IS  NOT  FLAT,  COMPUTE  RANGE  STEP 

c  •••  xt  MAY  BE  NECESSARY  TO  REDUCE  DZ  SUCH  THAT  DR  IS  SMALL  ENOUGH 

C  •••  TO  PRODUCE  ACCURATE  RESULTS 

IF ( ITYPEB.NE.3. AND. THETA. NE. 0.0)  DR  =  A  BSC DZ/ TAN( THETA ) ) 

C  •••  COMPUTE  DEPTH  INCREMENT  FOR  ADJUSTING  LAYERS 

DZZ=DR*TAN (THETA) 

C 

C  •••  GET  STARTING  FIELD 

IFC ISF.EC.O)  CALL  SF IELDC FRQ, CO , ZS . N . DZ . U ) 

IFCISF.NE.O)  CALL  UFIELD 
IFC IHNK.EQ.O)  GO  TO  69 
C 

C  •••  DIVIDE  STARTING  FIELD  BY  HANKEL  FUNCTION  IF  REQUESTED  BY  USER 

HNKsHNKLC XKO»RA) 

DO  63  1=1 ,N 
U(I)=U(I) / HNK 

63  CONTINUE 

69  CONTINUE 

C 

C  •••  COMPUTE  DEPTH  WRITE  INCREMENT  TO  NEAREST  DZ 

IF(WDZ.LT.DZ)WDZ=DZ 
IWZ=WDZ/DZ*. 5 
WDZ s IW  Z  *  DZ 

C 

C  •••  COMPUTE  DEPTH  PRINT  INCREMENT  TO  NEAREST  DZ 

IPZ=PDZ/DZ*.5 

IFC  PDZ.GT.O.O. AND. IPZ .EQ.O)  IPZ*  1 
PDZ=IPZ*DZ 
C 
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c 

72 

73 
7  4 

75 

76 

77 
80 


81 


<3*! 
3  5 
0 
C 


PRINT  PROBLEM  PARAMETERS 
WRITE (NPU, 72) 

F0RMAT(//1X,’IFD  SOLUTION') 

IFCISF.EQ.O)  WR ITE( NPU , 73 ) 

FORMAT( IX, 'GAUSSIAN  STARTING  FIELD') 

IF(ISF.NE.O)  WRITEC NPU, 74) 

FOHMAT( IX. 'USER  STARTING  FIELD') 

IF(IHNK.NE.O)  WRITEC  NPU, 75) 

FORMATC IX, ' STARTING  FIELD  DIVIDED  BY  HANKEL  FUNCTION') 
I F ( ITYPES. NE.O)  WR ITE ( N PU , 76 ) 

FORMATC IX, 'USER  SURFACE  CONDITION’) 

IFC ITYPEB.EQ. 1 )  WR ITE C N PU . 77 ) 

FORMATC IX. • USER  BOTTOM  CONDITION') 

WRITEC  NPU, 80)  FRQ.ZS.CO.RA.ZA.N 


FORMATC IX, 

'FRO 

s  ' , F7 . 1 '  HZ',/ 

C'ZS 

2 

'  ,F7. 

1 '  M' ,/ , IX, 

C'  CO 

2 

'  ,F7. 

1,'  M/SEC'./.IX, 

C’RA 

2 

'  ,F7. 

1 , '  M' ./ . IX, 

C'ZA 

= 

'  ,F7. 

1 , '  M' ,/ , IX. 

C'  N 

2 

•  .15) 

WRITEC  N PU , 8 1 )  DR.WDR 


FORMATC IX. 

•DR  3  ' ,F7 . 1 , ' 

C'WDR 

2 

' ,F7. 1 . '  M« ./ , IX. 

C  '  PDR 

2 

' ,F7. 1 . '  M' ./ . IX, 

C'DZ 

2 

' ,F7. 1 . '  M' ./ . IX. 

C’WDZ 

2 

' ,F7. 1 .'  M* ./ . IX, 

C*  PDZ 

2 

' ,F7. 1.'  M ' , / , 1 X , 

C ' ITYPEB 

2 

',1 i./.ix. 

C' ITYPES 

2 

’  . 11 ./ . IX, 

C ' RM AX 

_  • 

.  F8  .  1  ,  '  M'  ,/,/, 

C2X, ' LAYER 

MAX  DEPTHCM) 

M*  ,/.1X, 


DENSITY CG/CM* *3)  ATTC DB/WL) '  ./ ) 


DO  85  ILYR  sl.NLYR 
WRITEC NPU. 84) ILYR.ZLYRC ILYR) , RHOC ILYR) .BETAC ILYR) 
FORMATC 2X. I  3, 3 X, El 3.6,5X,E13.6,5X,E13.6) 

CONTINUE 

•••  PRINT  BOTTOM  DEPTHS  IF  REOUESTED 
IFC IBOT.EQ.O)  GO  TO  86 
THaTHETA"DEG 

WRITECNPU, 123)  R1.Z1.R2.Z2.TH 
IFC ISFLD. EO. 0 )  GO  TO  92 


•••  PRINT  STARTING  FIELD 
WRITEC  NPU. 87) 

87  FORMATC/ , IX, ' STARTING  FIELD') 

DO  90  I s 1 , N 
ZIsI'DZ 

WRITECNPU, 89)  I.ZI.UCI) 

89  FORMATC 1X.I4,3X,F10.2,3X,'( ' ,E12.5,2X,E12.8,'  )’) 

90  CONTINUE 

92  CONTINUE 

C 

IFC ISVP.EO.O)  GO  TO  96 
C  •••  PRINT  SVP 
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WRITE(NPU,93)  RA 

93  F0RMAT(//1X, ' SOUND  VELOCITY  PROFILE  AT  RANGE  '.F8.1.'  METERS',/) 

DO  95  Isl.NSVP 

WRITE  (NPU, 9*0  I.ZSVP(I) .CSVP(I) 

9«  FORMAT(1X.Itt,3X,F8.1,3X,F8.1) 

95  CONTINUE 

96  CONTINUE 
C 

C  *•*  COMPUTE  'B'  COEFFICIENT 

BCOFsCMPLX(0.0, .5/XKO) 

C 

C  •••  ASSIGN  OUTPUT  FILE  IF  OUTPUT  REQUESTED 

IF(WDR.EQ.O.O)GO  TO  98 
CALL  ASSIGN(NOU.'IFD.OUT') 

c 

C  WRITE  SELECTED  PARAMETERS  FOLLOWED  BY  STARTING  FIELD 

WRITE(NOU)FRQ,ZS,CO,ISF.RA,ZA,N,IHNK,ITYPEB,ITYPES.RMAX,DR,WDR,DZ. 

CNLYR.ZLYR.RHO.BETA 

MZsN/IWZ 

WRITE(NOU)  N  Z , R A , WDZ , ( U ( I ) . IsIWZ.N.IWZ) 

C 

98  CONTINUE 


c 

III 

XPR  s 

INITIALIZE 

RA+PDR 

RANGE 

VARIABLE 

AT 

WHICH 

SOLUTION 

IS 

TO 

BE 

PRINTED 

c 

»•* 

INITIALIZE 

RANGE 

VARIABLE 

AT 

WHICH 

SOLUTION 

IS 

TO 

BE 

RECORDED 

XWRsRA+WDR 

IF(WDR.EQ.O.O)  XWR=RA*RMAX*1 .0 
C 

C  •••  SAVE  RANGE  STEP 

OLDRsDR 

C 

C  INITIALIZE  PARAMS  FOR  DIAG  THEN  COMPUTE  MAIN  DIAGONALS  X  AND  Y 

Nisi 
DRIsDR 

C  COMPUTE  X  DIAGONAL 

CALL  DIAG 

C  •••  COMPUTE  Y  DIAGONAL 

CALL  DIAG 

C 

C  MAIN  LOOP  STARTS  HERE  ••* 

C  SOLUTION  WILL  BE  ADVANCED  FROM  RANGE  RA  TO  RANGE  RA*DR 

C 

100  CONTINUE 
IDIAGsO 
C 

C  •••  DOES  SVP  CHANGE  BEFORE  BOTTOM  PR*. ILE? 

IF( RSVP.GE. R2)  GO  TO  105 

C  •••  YES.  DOES  SVP  CHANGE  BEFORE  NEXT  SOLUTION  RANGE? 

:f(ra*df.le.rsvp)  go  to  131 

C  **•*  YES.  ADJUST  DR  SO  THAT  SOLUTION  ADVANCES  TO  RSVP 

DRsRSVP-RA 
GO  TO  126 

C  •••  BOTTOM  PROFILE  CHANGES  BEFORE  OR  AT  SAME  RANGE  AS  SVP 
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105  CONTINUE 

•••  DETERMINE  IF  BOTTOM  DEPTHS  ARE  TO  BE  UPDATED 
•••  FIRST  PASS  -  RA  WILL  BE  .LT.  R2 
•••  UPDATE  BOTTOM  DEPTHS? 

IF( RA.GE. R2)  GO  TO  110 

•••  NO. 

•••  DETERMINE  IF  RANGE  STEP  TOO  LARGE 
IF( RA+DR . LE. R2)  GO  TO  131 

•••  RANGE  STEP  IS  TOO  LARGE  -  RESET  DR  -  ADV  SOLUTION  TO  R2 

DR=R2-RA  !  . 

GO  TO  126 

•••  UPDATE  BOTTOM  DEPTHS 
10  CONTINUE 
R  1  sR2 
Z 1  sZ2 

ITRK*ITRK+1 
R2aTRACK( ITRK. 1 ) 

Z2=TRACK( ITRK.2) 

•••  TWO  DEPTHS  AT  SAME  RANGE  INDICATE  VERTICAL  DISCONTINUITY. 

ADVANCE  TRACK  FORWARD . 

IF(R1.EQ.R2)  GO  TO  110 

•••  RESTORE  DR 
DRzOLDR 

COMPUTE  SLOPE  OF  BOTTOM 
THETAsATAN2( Z2-Z 1 .R2-R1) 

IF  BOTTOM  IS  NOT  FLAT.  COMPUTE  NEW  RANGE  STEP 
IF(THETA.EQ.O)  GO  TO  120 

•••  IF  BOTTOM  OF  ARTIFICIAL  LAYER  IS  FLAT,  DO  NOT  RECOMPUTE  DR. 
IF( ITYPEB.EQ. 3)  GO  TO  120 
DR=ABS(DZ»COS< THETA )/SIN< THETA) ) 

•••  IF  RANGE  STEP  TOO  LARGE,  PRINT  WARNING  MESSAGE.  RECOMPUTE  DR. 
IF( RA+DR. LE. R2)  GO  TO  120 
W  RITE ( N  PU , 1 15) 

115  FORMATt IX, ' RANGE  STEP  TOO  LARGE  FOR  BOTTOM  IRREGULARITIES') 

C  STOP  ! . 

DR.R2-RA  !  . 

C 

120  CONTINUE 
C 

C  •••  PRINT  BOTTOM  DEPTHS 

IF( I30T.E0.0)  GO  TO  126 
THsTHETA*  DEG 

WRITE(NPU. 123)  R1.Z1.R2.Z2.TH 
123  F0RMAT(//1X, 'BOTTOM  DEPTHS  ',//,1X, 

C'RI  s  * , F8 . 1 , '  M ' , / , 1 X , 

C'Zl  1  '  ,F8. 1  ,  '  M'  ,/ , 1 X, 


A-14 


TR  6659 


i 


C '  R  2  *  ' , F8  .  1  .  '  M'  ./  , 1 X, 

C'Z2  *  ' ,F8. 1 . '  M*  ,/ , IX, 

C ' THETA  *  ’ ,F8. 1 . ’  DEG' ) 

126  CONTINUE 

DZZsDR»TAN(THETA) 

WRITE(NPU, 128)  DR , R A 

128  FORMAT( IX, ' RANGE  STEP  s  '.FU.O,'  M  AT  RANGE  ' , F8 . 1 , '  M' ) 

C 

C  •••  ADVANCE  RANGE  ONE  STEP 

130  RAsRA*DR 
IDIAGs 1 
GO  TO  132 

131  CONTINUE 
RAsRA+DR 

132  CONTINUE 

C  READ  NEW  SVP  PROFILE  FLAG? 

IF(RA.GE.RSVP)  R EA D( N IU , • , EN D= 3 3 >  KSVP 
C 

C  •**  IF  KSVP  NOT  0,  SUBROUTINE  USVP  IS  CALLED  FOR  SVP.  USER  IS 

C  •••  RESPONSIBLE  FOR  ADJUSTING  DEPTHS  OF  LAYERS  WHEN  BOTTOM 

C  IS  NOT  FLAT. 

IF(KSVP.EQ.O)  GO  TO  13“ 

CALL  USVP 
IDIAGs 1 
GO  TO  148 
C 

C  •••  NEW  SVP  AT  ADVANCED  RANGE? 

134  IF(RA.GE.RSVP)  GO  TO  145  ! . 

C 

C  •••  NO.  IS  BOTTOM  FLAT? 

135  IF(THETA.EO.O.O)  GO  TO  166 
MLYRsNLYR 

IF( ITYPEB.EQ. 3)  MLYRsNLYR-1 
C 

C  •••  UPDATE  DEPTHS  OF  LAYERS 

C  •••  ASSUMES  LAYERS  ARE  PARALLEL  AND  FOLLOW  BOTTOM  CONTOUR 

C  IF  ITYPEB  s  3,  BOTTOM  OF  ARTIFICIAL  LAYER  REMAINS  FLAT. 

DO  138  ILYR* 1 , MLYR 
ZLYR( ILYR) sZLYR< ILYR)*DZZ 

138  CONTINUE 

C  N1*ZLYR(  D/DZ 

Nisi 

ZAsZLYR(NLYR) 

IF( ITYPEB.EQ. 3. AND. ZLYR(NLYR) .LT.ZLYR(  MLYR) ) WRITE (NPU,  1  ?9) 

139  FORMAT ( 1 X , ' ER R .  DEPTH  OF  ARTIFICIAL  LAYER  LESS  THAN  LATER  ABOVE’) 
C 

C  ADJUST  DEPTHS  OF  PROFILES  IN  SLOPING  LAYERS 

C  •••  ASSUMES  SAME  SVP  IN  LAYERS 

NWSVPs I XS VP( 1 ) 

DO  140  IsNWSVP* 1 , NSVP 
ZSVP(I)sZSVP(I)*DZZ 

140  CONTINUE 
GO  TO  167 

C 
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C  •••  GET  NEXT  SVP 

145  CONTINUE 

CALL  SVP(NLYR.ZLYR.RHO,BETA,IXSVP.NSVP,ZSVP,CSVP) 

IDIAGs  1 

148  CONTINUE 

C  ERROR  DETECTED? 

IF( N  SVP . EQ. 0 )  GO  TO  33 
IF( ITYPEB. NE . 3)  ZAsZA+DZZ 
C  •••  NO.  READ  RANGE  OF  NEXT  PROFILE? 

IF( RA.GE. RSVP)  R E AD( N IU . • . EN Ds 1 50 )  RSVP 

149  CONTINUE 

C  •••  ARTIFICIAL  ABSORBING  LAYER? 

IF( ITYPEB.LT. 2)  GO  TO  155 
IF( ITYPEB. GT. 3)  GO  TO  155 
C  •••  YES.  UPDATE  DENSITY.  ATTEN  AND  SPEED. 

IF( ZA.LT.ZLYR(NLYR) )  GO  TO  155 

NLYRsNLYR* 1 

ZLYRC  NLYR) sZ A 

RHO( NLYR) sRHO< NLYR- 1 ) 

BETA(  NLYR) sB ET A<  N LYR- 1 ) 

NSVPsNSVP* 1 
IXSVP(NLYR)sNSVP 
ZSVP(NSVP)sZ LYR (NLYR) 

CSVP(NSVP) sCSVP( NSVP-1 ) 

GO  TO  155 

C  •••  SET  RSVP  LARGE  SO  THAT  LAST  PROFILE  IS  USED  FOR  REMAINDER  OF  PROBLEM 

150  RSVP=1.0E*38 
GO  TO  149 

C 

C  PRINT  SVP 

155  IF( ISVP.EO.O)  GO  TO  165 
WRITE (NPU, 93)  RA 
DO  160  Is  1 , NSVP 

WRITE( NPU . 9*0  I.ZSVP(I) .CSVP(I) 

160  CONTINUE 

165  CONTINUE 
C 

C  IF  NEW  DR  AND/OR  SVP  -  UPDATE  X  AND  Y  DIAGONALS 

166  IF( IDIAG.EQ.O)  GO  TO  170 
Nisi 

167  CALL  DIAG 
DRIsDR 

170  CONTINUE 
C 

C  •••  ADVANCE  SOLUTION  ONE  STEP  FORWARD 

NOLDsN 
CALL  CRNK 
C 

C  APPLY  ABSORPTION  IF  ITYPEB  i  2  OR  i 

IF( ITYPEB.LT. 2)  GO  TO  185 
IF( ITYPEB. GT. 3)  GO  TO  185 
MMsZLYRC  NLYR-1 ) /DZ 
NAsN-MK 

IF(NA.GT.O)  GO  TO  175 
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IF( N A . EQ . 0 )  GO  TO  185 
WRITE!NPU, 174)  RA 

174  FORMAT! 1 X, • ERR  IN  THICKNESS  OF  ABSORBING  LAYER  AT  ',F8.1,'  M') 
STOP 

175  CONTINUE 

C  SEE  AESD  PE  MODEL  BY  BROCK  -  NORDA  TECH  NOTE  12  -  JAN  78 

DO  180  1=1. NA 

ATT=EXP(-.01»DR»EXP(-(<I-NA)/(NA/3.0))#,2.0)) 
U(MM+I)=U!MM+I)*ATT 
180  CONTINUE 
185  CONTINUE 
C 

C  TERMINATE  RUN  IF  N  HAS  BEEN  DECREMENTED  TO  5  -  ARBITRARY 

IF( N.GT.5)  GO  TO  200 
W  RITE ( N  PU , 190) 

190  FORMAT! IX, ' PROGRAM  TERMINATED  -  ONLY  FIVE  POINTS  REMAIN') 

GO  TO  400 
200  CONTINUE 
C 

C  •••  TERMINATE  RUN  IF  N  HAS  BEEN  INCREMENTED  BEYOND  MAXIMUM 

IF(N.LE.MXN)  GO  TO  250 
WRITE(NPU.210)RA 

210  FORMAT! IX. 'MAX  N  EXCEEDED  AT  RANGE  '.F10.2.'  M.') 

GO  TO  400 
250  CONTINUE 
C 

C  •••  DETERMINE  IF  NEW  POINT  ADDED 

IF! N . LE . SOLD )  GO  TC  255 
NIsNOLD 
CALL  DIAG 
DR  1 =DR 

255  CONTINUE 
C 

C  •••  IF  SOLUTION  IS  TO  BE  WRITTEN  ON  DISK, 

C  WRITE  SELECTED  PRESSURE  FIELD  AT  RANGE  RA 

IF(XWR.GT.RA)  GO  TO  260 

WRITE! NOU)  HZ , R A , WDZ , ! U! I) , I  * IWZ , N , IWZ) 

C 

C  •••  DETERMINE  NEXT  RANGE  AT  WHICH  TO  WRITE  SOLUTION  ON  DISK 

XWRs  XWR*WDR 
260  CONTINUE 
C 

C  •••  DETERMINE  IF  SOLUTION  IS  TO  BE  PRINTED 

IF! XPR.GT. RA.OR . IPZ . EQ. 0. OR. PDR .EQ. 0.0)  GO  TO  ?50 
C 

C  •••  PRINT  RANGE 

RTEMP=  R A/  1  000 . 0 
wr:te:npu.270)  rtemp 

270  FORMAT! /5X, ' RANGE  =',E15.8,'  KM.  •/) 

C 

C  •••  COMPUTE  HANKEL  FUNCTION 

IF! IHNK.NE.C)  HNKsHNKL! XKO*RA) 

C 

C  •••  COMPUTE  AND  PRINT  PROPAGATION  LOSS  AT  EACH  IPZ'TH  DEPTH 
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WRITE! NPU, 275) 

275  FORMAT! 6 X, ' I • , 9X , • Z( I ) • . 6X , • LOSS( DB) ’  .  1  4X . • 0 ( I)  • ) 

DO  300  IsIPZ.N.IPZ 
Z Is  I *DZ 
PLsCABS ( U ( I ) ) 

IF! IHNK.EQ. 1)  PLsCABS!U!I)»HNK) 

IF( PL.LE.0.0)  GO  TO  288 
PLs-20.0*A LOGIC! PL) 

I F( IHNK.EQ.O)  PLsPL*10.0»ALOG10!RA> 

GO  TO  289 

288  PLs999.9 

289  WRITE( NPU , 295)  I.ZI.PL.U(I) 

295  FORMAT!  2X, 15.  (  SX.FIO^^X.FIO^J.’X.'C.EIP.S^X.E^.S.’  )  ’  ) 

300  CONTINUE 
C 

C  •••  DETERMINE  NEXT  RANGE  AT  WHICH  TO  PRINT  SOLUTION 

XPRs  XPR  +  PDR 

C 

C  •••  TERMINATE  RUN  IF  SOLUTION  AT  MAXIMUM  RANGE  HAS  BEEN  OBTAINED 

350  IF( RA.LT. RMAX)  GO  TO  100 

C 

C  •••  TERMINATE  run 

400  IF(WDR.NE.O)CALL  CLOSE(NOU) 

W  R ITE( NPU, 10)  TM 

CALL  TIME! TM) 

WR ITE ! N  PU , 10)  TM 
WR ITE ! N  PU , 40 1 ) 

401  FORMAT! IX. ' END  OF  RUN') 

STOP 

END 
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SUBROUTINE  DIAG 


*  COMPUTES  RANGE  AND  DEPTH  DEPENDENT  MAIN  DIAGONALS  OF  MATRICES  * 

*  DIAG  IS  CALLED  WHENEVER  BOTTOM  DEPTH  OR  SVP  CHANGES  * 


***  SEE  MAIN  PROGRAM  FOR  DEFINITIONS 
***  SUBROUTINE  DIAG  RETURNS: 

ACOFX  -  COEFFICIENT  'A'  AT  BOTTOM  -  AT  RANGE  RA 
ACOFY  -  COEFFICIENT  *A'  AT  BOTTOM  -  AT  RANGE  RA-DR 
BTA  -  ARRAY  -  PARTIAL  SOLUTION  OF  SYSTEM  OF  EQUATIONS 
X  -  ARRAY  -  MAIN  DIAGONAL  OF  MATRIX  AT  RANGE  RA 

Y  -  ARRAY  -  MAIN  DIAGONAL  OF  MATRIX  AT  RANGE  RA-DR 


PARAMETER  MXLYR-10 1 , MXN-10000 , MXSVP-101 , MXTRK-1 0 1 ,NIU«1 , 

C  NOU*2,NPU«8 

COMPLEX  ACOFX , ACOFY , BCOF , BOTX , BOTY , BTA , HNK , HNKL , SURX , SURY , TEMP , 

C  U  ,  X ,  Y 

COMMON  /TFDCOM/ACOFX, ACOFY, ALPHA, BCOF, BETA(MXLYR) , BOTX, BOTY, 

C  BTA(MXN) , CO ,CSVP (MXSVP) , DR , DR1 , DZ , FRQ , t HNK , ISF , ITYPEB , 

C  I TYPES , IXSVP ( MXLYR) , KS VP , N , N 1 , NLYR , NSVP , NWS VP , R 1 2 ( MXN ) , RA , 

C  RHO ( MXLYR) , RSVP, SURX, SURY, THETA, TRACK (MXTRK, 2) ,U(MXN) , 

C  X ( MXN)  , XKO , Y  ( MXN )  , ZA , ZLYR (MXLYR)  , ZP , ZS , ZSVP (MXSVP) 

DATA  PI/3. 14 1592554/, DEG/57. 29578/ 

COMPLEX  B , XN 1 , XN 2 , EYE 
DATA  IEX/0/ 

***  COMPUTE  NEW  X  AND  Y  DIAGONALS 
EYE=CMPLX  (0 . 0 , 1 . 0 ) 

***  GET  INDICES,  DENSITY,  AND  ATTENUATION  FOR  FIRST  LAYER 

ILYR-1 

L»1 

M=»IXSVP(1) 

Rl=RHO(l) 

BETA1 *BETA ( 1 ) 

DO  10  I *N 1 ,  N 

***  TRANSFORM  X  INTO  Y 

Y(I)  *-X ( I ) -EYE *2 .0*DZ*DZ*XKO* (1 . 0  +  R 1 2 (I) ) * (1 . 0/DR1+1 .O/DR) 

0  CONTINUE 

DO  100  I-N1,N 

***  ZI  IS  RECEIVER  DEPTH 

ZI»I*DZ 

***  IS  RECEIVER  IN  THIS  LAYER? 

0  IF  ( Z I .  LE.ZLYR  (ILYR)  )  GO  TO  <19 

***  NO.  ''LL  LAYERS  CHECKED? 

IF(ILYR.SO.NLYR)  GO  TO  52 

C  ***  NO.  SET  INDICES,  DENSITY,  AND  ATTENUATION  FOR  NEXT  LAYER. 

ILYR-ILYR+1 
L-M+l 

M« IXSVP (ILYR) 

Rl*RHO (ILYR) 

B ETA 1=B ETA (ILYR) 
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GO  TO  20 

C  ***  OEPTH  ZI  IS  IN  THIS  LAYER. 

49  CONTINUE 

C  ***  DETERMINE  SOUND  SPEED  Cl  AT  DEPTH  ZI 

DO  50  J  =  L,M 

IF ( Z I .GT . ZSVP ( J ) )  GO  TO  50 
L=J 

C  ***  INTERPOLATE 

CI=CSVP (J-l ) + (CSVP ( J ) -CSVP ( J-l ) ) * ( Z I-ZSVP ( J-l ) ) / (ZSVP (J ) -ZSVP (J-l ) 
C) 

GO  TO  <50 

50  CONTINUE 

C  ***  EXTRAPOLATE 

52  CONTINUE 

I F ( ZSVP (M) . NE . ZSVP ( M- 1 ) ) GO  TO  53 
CI»CSVP(M) 

GO  TO  54 

53  CI=CSVP(M-1)+ (CSVP(M) -CSVP(M-l) ) * ( Z I-Z SVP (M-l ) )/ 

C ( ZSVP ( M) -ZSVP ( M-l ) ) 

54  IF ( I  EX. EQ. 1 )  GO  TO  50 
WRITE (NPU , 56 ) 

55  FORMAT (IX, 'WARNING.  EXTRAPOLATION  OF  SVP  PERFORMED.') 

I  EX3 1 

50  CONTINUE 

C  ***  SAVE  SPEED  IN  MEDIUM  1 

C1=»CI 

C  ***  IS  RECEIVER  ON  OR  WITHIN  1  DZ  OF  INTERFACE? 

I F ( Z LYR ( I LYR) -ZI.GE.DZ)  GO  TO  55 
C  ***  YES.  IS  THIS  BOTTOM  LAYER? 

I F ( I LYR . EQ.NLYR)  GO  TO  55 
C  ***  NO.  GET  PARAMETERS  FOR  MEDIUM  2 

R2=RHO ( I LYR+1 ) 

BETA2=*BETA  (  T  LYR+1 ) 

C2 =»CSVP ( I XSVP ( I LYR ) + 1 ) 

GO  TO  ?0 
55  CONTINUE 

C  ***  NOT  AN  INTERFACE.  USE  MEDIUM  1  PARAMETERS. 

C2*C1 

R2*R1 

BETA?*BETA1 

C  ***  COMPUTE  DENSITY  RATIO 

70  R 12 ( I ) *R 1/R? 

C  **•  THE  NEXT  FEW  LINES  COMPUTE  THE  X  OIACONAL 

XN-C0/C1 

I F ( BETA1 . LT . 0 . 0)  BETA1 »ALPHA*C1/FR0 

XN1-CMPLX (XN*XN,XN*XN*BETAt/27.2B75?7) 

XN*C0/C7 

IF (PETA2.LT. 0.0)  BETA2*ALPHA*C7/FRO 

XN2*CMPLX  (XN*XN,XN*XN*BETA2/?7.2R75?-7) 

B-  (-XK0* . 5 * (  (XN1-1 . ) +R12 (I) * (XN2-1 .0) ) ) 

X (I) »DZ*DZ*XX0« (B-CMPLX(0.0,  +  2.0* (1.0+R12 (I)  )  /DR)  ) +1+R12 (I) 

100  CONTINUE 

C 

C  ***  COMPUTE  BTACI1)  THROUGH  BTA(N) 

C  ***  ARRAY  BTA  CONTAINS  PARTIAL  SOLUTION  of  SYSTEM  OF  EQUATIONS 

IF(Nl.EO.l)  GO  TO  105 
«*N1 

GO  TO  1^5 
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105  BTA ( 1) “X (1 ) 

M-2 

105  DO  110  I=M,N 

BTA(I)»X(I)-R12(I-1)/BTA(I-1) 

10  CONTINUE 

***  COMPUTE  'A'  COEFFICIENTS  AT  BOTTOM 
ACOFY-ACOFX 

ACOFX-CMPLX (0 . 0 ,0 . 5) *XK0* (XN2-1 .0) 

RETURN 

END 
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SUBROUTINE  CRNK 

C  IKKKKKKKKKKK KIKKKIKI 

C  MiMMIIIMtMItltXKMIMMMIMIIIMtlKItliMMlilflMMtlfff 

C  •  THIS  ROUTINE  USES  THE  CR ANK-N ICOLSON  METHOD  FOR  SOLVING  THE  * 

C  •  SYSTEM  OF  EQUATIONS.  THE  SOLUTION  IS  ADVANCED  FROM  RANGE  RA-DR  • 

C  •  TO  RANGE  RA  WHERE  RA  IS  THE  ADVANCED  RANGE.  • 

C  KKKKKKKKK KKKKKKKIKKKKKKK  KKKKKKKKK | 

C  IKKKKKKIKIKKKKKKKKKKIIIKKIIIKKKIKIKIIKKI 

C  •••  SYSTEM  OF  EQUATIONS  IS  AS  FOLLOWS: 

C 


c 

ALL  VALUES 

ARE  AT  ADVANCED  RANGE 

c 

III 

in 

III  IK 

in  ••• 

c 

*  X(1) 

-R 1 2(  1  ) 

0  » 

•  U(  1  )  • 

•  SURX  * 

c 

•  -1 

X  (  2 ) 

-R 1 2 ( 2)  •  . 

•  U( 2)  •  - 

•  0  • 

c 

•  0 

-1 

X(N)  • 

•  U  (  N )  • 

•  BOTX  • 

c 

IK 

••• 

••• 

c 

c 

ALL  VALUES 

ARE  AT  PRESENT 

1  RANGE 

c 

III 

in  ••• 

c 

•  Yd) 

♦  R 1 2(  1  ) 

0  • 

•  ud>  • 

•  SURY  • 

c 

•  *1 

Y  (  2  > 

♦R 1 2 ( 2)  •  . 

*  U  (  2)  •  ♦ 

•  0  • 

c 

•  0 

4-1 

Y(N)  • 

•  U  (  N )  • 

•  BOTY  • 

c 

III 

in 

m  ••• 

C 


C  D  ARRAY  -  RIGHT  HAND  SIDE 

C  RA  RANGE  TO  WHICH  SOLUTION  IS  TO  BE  ADVANCED  -  METERS 

C  -  RA  IS  INCREMENTED  BY  DR  PRIOR  TO  ENTERING  THIS  ROUTINE 

C  TA  TEMPORARY  VARIABLE  USED  IN  MANIPULATION  OF  BOTTOM  POINT 

C  TB  TEMPORARY  VARIABLE  USED  IN  MANIPULATION  OF  BOTTOM  POINT 

C 

C  •••  SEE  MAIN  PROGRAM  FOR  OTHER  DEFINITIONS 

C  •••  UNDEFINED  VARIABLES  ARE  TEMPORARY  VARIABLES 

C  •••  SUBROUTINE  CRNK  RETURNS: 

C  N  NUMBER  OF  EQUI-SPACED  POINTS  IN  U  AT  RANGE  RA 

C  U  ARRAY  -  COMPLEX  ACOUSTIC  PRESSURE  FIELD  AT  RANGE  RA 

C  -  INCLUDES  BOTTOM  POINT  -  DOES  NOT  INCLUDE  SURFACE  POINT 

C  •••  SUBROUTINE  TRID  IS  CALLED  TO  SOLVE  THE  TRIDIAGONAL  MATRIX 

C  KKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK 

C 

PARAMETER  MXLYRslOI ,  MXNs  1 0000 , MXSVPr 1 0  1 . MXTRKslOI , NIUsI  . 

C  NOU*2,NPUs6 

COMPLEX  ACOFX. ACOF Y , BCOF , BOTX , BOTY , BTA , HNK , HNKL . SURX , SUR Y , TEMP . 

C  U.X.Y 

COMMON  /IFDCOM/ACOFX, ACOFY, ALPHA. BCOF, BETA(MXLYR) .POTX.BOTY, 
BTA(MXN) , CO.CSVP(MXSVP) . DR.DR1 , DZ.FR Q.IHNK.TSF, TTYPEB , 

I  TYPES, IXSVP(MXLYR) .KSVP, N,N1 . NLYR , NS VP . NWS VP . R 1 2 ( MXN )  ,  RA. 
RHO(MXLYR) . RS V P , SUR X . SUR Y . TH ET A , TR ACK( MXTR K . 2 ) ,  U(MXN) . 

X(MXN) .XKO.Y(MXN) , ZA , ZLYRC  MXLYR) , ZP , ZS , ZS VP ( MXS VP ) 

DATA  PI/ 3. 1« 1592654/ .DEG/57.29578/ 

DIMENSION  D(MXN) 

COMPLEX  D, SUM , A , EX . P , 0 . XX . S 1 . S2 , C A . CB . CC , CD . TA , TB 
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c 

C  •••  CALL  SCON  FOR  SURFACE  CONDITION 

CALL  SCON 

»•*  COMPUTE  RIGHT  HAND  SIDE  DM)  THROUGH  D(N-1) 

•••  D(N)  IS  COMPUTED  LATER 

D< 1 ) sY( 1 )»U( 1 )+R 12( 1 ) »U(2)+SURY*SURX 

DO  5  Is2 , N- 1 

D( I)sU( I-1)+Y(I)»U(I)+R12( I)*U(I+1) 

BOTTOM  TYPE 
IBs ITY  PEB+  1 

GO  TO  ( 10,40,80,82)  .IB 

BOTTOM  IS  RIGID  -  ITYPEB  s  0 
0  CONTINUE 

IF( THETA . GT. 0 . 0 )  GO  TO  30 
IF ( THETA . NE . 0 . 0 )  GO  TO  15 

•••  RIGID  BOTTOM  IS  FLAT 
BOTYsU ( N ) 

D(N) sU(N-1)*Y(N)»U(N)*BOTY 
TAsO.O 

•••  BOTXsU ( N )  AT  ADVANCED  RANGE 
TBs  1 . 0 

CALL  TRID(X.U,D,BTA.R12,N,N,TA,TB) 

RETURN 

•••  RIGID  BOTTOM  SLOPES  UPWARD  -  DELETE  A  POINT  -  USE  2ND  ORDER  ODE 
15  BOTYsU(N) 

NsN-1 

D(N)sU(N-1)-*Y'-<)*U(N)*BOTY 
COTsCOS( THETA )/ SIN (THETA) 

Ps-COT/BCOF 

Qs(ACOFX*CMPLX(O.O.XKO))/BCOF 
XX  sCSQRT ( P*  P-4 . 0*0) 

S1*(-P*XX)/2.0 
S2s ( -P-XX) / 2 . 0 
CAs ( 1.0-S1*DZ)/(-XX»DZ) 

CBs1.0/(XX»DZ) 

CCs 1 . 0-C A 
CDsCB 

TAs-CD»CEXP(S1»DZ) ♦CB*CEXP(  S2*DZ) 

TBsCC»CEXP(S1»DZ) *CA»CEXP( S2»DZ) 

CALL  TRID(X,U,D,BTA,R12.N,N,TA,TB) 

RETURN 

C 

30  CONTINUE 

C  •••  RIGID  BOTTOM  SLOPES  DOWNWARD  -  USE  2ND  ORDER  O.D.E. 

COTsCOS( THETA) /SIN ( THETA) 

Ps-COT/BCOF 

Qs(ACOFY+CMPLX(O.OtXKO) )/BC0F 
XXsCSORTC  P*P-4.0*Q) 

S1s(-P*XX)/2.0 
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S2s ( -P-XX) /  2 . 0 

CAs( 1 .0-S1»DZ)/(-XX»DZ) 

CBs 1 .0/ (XX»DZ) 

CCsl .0-CA 
CDsCB 

TAs-CD*CEXP(S1#DZ)+CB#CEXP(S2*DZ) 

TBsCC»CEXP(S1»DZ) ♦CA^CEXPC S2»DZ) 

MODIFY  LOWER  AND  Y  DIAGONALS  IN  ROW  N  BY  TA  AND  TB  -  THEN 
•••  COMPUTE  D<  N ) 

D(N)*( 1 .0*TA)»U(N-1)*(Y(N)+TB)«U(N) 

•••  COMPUTE  TA  AND  TB  FOR  LOWER  AND  X  DIAGONALS  IN  ROW  N 
Qs  ( ACOFX  +  CM  PLX ( 0.0, XKO) )/BC0F 
XXsCSQRT( P»P-4.0»Q) 

S1s(-P*XX)/2.0 
S2s(-P-XX)/2.0 
CAs( 1 .0-S1«DZ)/(-XX«DZ) 

CBs 1 . 0/  ( XX*DZ) 

CCs 1 .0-CA 
CDsCB 

TAs-CD*CEXP(S1*DZ)+CB,CEXP(S2#DZ) 

TBsCC»CEXP(Sl»DZ) *C A«CEXP( S2»DZ) 

CALL  TRID(X,U.D.BTA,R12,N,N.TA,TB) 

*»*  ADD  A  POINT 

CBs< (U(N)-U(N-I) )/DZ-S1»U(N) )/(-XX) 

CAsU(N)-CB 

NsN+1 

U( M) sCA*CEXP(S  1»DZ) ♦CB*CEXP(S2»DZ) 

RETURN 

CONTINUE 

BOTTOM  CONDITION  SUPPLIED  BY  USER  -  ITYPEB  s  1 
IF( THETA. GT. 0.0)  GO  TO  60 
IF( THETA . NE . 0 . 0)  GO  TO  70 

FLAT  BOTTOM 

•••  USER  SUPPLIES  BOTY  ON  BOTTOM  INTERFACE  AT  PRESENT  RANGE 

•••  USER  SUPPLIES  BOTX  ON  BOTTOM  INTERFACE  AT  ADVANCED  RANGE 

CALL  BCON 

TAsO.O 

TBsO. 0 

NsN-1 

D(N)sU(N-1)«.Y(N)»U(N)+B0TY*B0TX 
CALL  TRID(X,U,D,BTA,R12,N,N,TA,TB) 

NsNt-1 

U ( N ) sBOTX 

RETURN 


0 
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C  BOTY  -  ONE  DZ  IN  BOTTOM  AT  PRESENT  RANGE 

C  BOTX  -  ON  BOTTOM  INTERFACE  AT  ADVANCED  RANGE 

CALL  BCON 
TAsO.O 
TBsO.O 

D(N)sU(N-1)+Y(N)»U(N)+BOTY*BOTX 
CALL  TRID(X,U,D.BTA.R12.N,N,TA,TB) 

Ns  N  ♦  1 

U ( N ) sBOTX 

RETURN 

C 

70  CONTINUE 

LAYER  SLOPES  UP 
•••  USER  SUPPLIES  : 

BOTX  -  ONE  DZ  BELOW  BOTTOM  INTERFACE  AT  ADVANCED  RANGE 
CALL  BCON 
NsN-1 

D(N)sU(N-1)*Y(N)*U(N)+BOTY*BOTX 
TAsO.O 
TBsO.O 

CALL  TRIDCX.U.D. BTA , R 1 2 . N , N , TA , TB) 

RETURN 


•••  I TYPES  s  2  —  ARTIFICIAL  ABSORBING  LAYER  INTRODUCED 
0  CONTINUE 

IF( THETA . NE . 0 . 0)  GO  TO  85 

•••  BOTTOM  OF  ABSORBING  LAYER  IS  FLAT 
•••  ITYPEB  s  2  OR  3  —  ARTIFICIAL  ABSORBING  LAYER 
2  U( N ) sO  .  0 

TAsO.O 
TBsO.O 

CALL  TRID(X.U.D.BTA.R12.N-1 .N.TA.TB) 

RETURN 

###  ITYPEB  s  2  —  ARTIFICIAL  ABSORBING  LAYER  NOT  FLAT 

5  IF( THETA . GT. 0 . 0 )  GO  TO  90 

•**  BOTTOM  OF  ABSORBING  LAYER  SLOPES  UP 

TAsO.O 

TBsO.O 

U(N-1)=0.0 

U(N)s0.0 

CALL  TRID(X,U.D,BTA,R12.N-2,N-1 ,TA,TB) 

DELETE  A  POINT 
NsN-1 
RETURN 
C 

90  CONTINUE 

C  •••  BOTTOM  OF  ABSORBING  LAYER  SLOPES  DOWN 

D(N) sU(N-1 ) 

TAsO.O 
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TBsO.O 

CALL  TRID(X.U,D,BTA,R12,N,N,TA,TB) 
C  »•«  ADD  A  POINT 

NsN  +  1 
U(N)=0.0 
RETURN 
END 
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SUBROUTINE  TR ID(  X , U , D , BTA , R 1 2 , L , M . TA , TB) 

C 

C  •••  SPECIALIZED  VERSION  OF  METHOD  PRESENTED  ON  PG  442  OF  CARNAHAN 

C  •••  ET  AL  FOR  SOLVING  A  TRIDIAGONAL  MATRIX. 

C  •••  TRID  RETURNS  THE  SOLUTION  FIELD  IN  U 

C  •••  LOWER  DIAGONAL  =  -1 

C  •••  UPPER  DIAGONAL  =  -R12 

C  BTA  IS  PARTIAL  SOLUTION  COMPUTED  IN  ROUTINE  DIAG 

C  •••  D  CONTAINS  R.H.S. 

C  •••  GAMMA  -  TEMPORARY  STORAGE 

C  •••  TA  -  TEMPORARY  VARIABLE  USED  IN  MANIPULATION  OF  BOTTOM  POINT 

C  »**  TB  -  TEMPORARY  VARIABLE  USED  IN  MANIPULATION  OF  BOTTOM  POINT 

C 

c 

PARAMETER  MXNslOOOO 

DIMENSION  X( 1 ) ,U( 1 ) , GAMM A( MXN5 . D( 1 ) , BTA( 1 ) , R12( 1 ) 

COMPLEX  X,U, GAMMA. D. BTA, TA.TB 

C  X  AND  LOWER  DIAGONAL  IN  ROW  L  MUST  BE  MODIFIED  BY  TB  AND  TA 

BTA(L)=(X(L)-TB)-( ( 1 . O+T A) • R 1 2 ( L- 1 > > / BT A( L- 1) 

GAMMA< 1)=D( 1 ) / BTA( 1) 

DO  10  1=2. L 

GAMMA(  I) =(D( I)+GAMMA( 1-1 ) )/BTA( I) 

10  CONTINUE 

GAMMA(  L) =GAMMA( L ) ♦ ( TA»G AMMA< L- 1 ) )/BTA( L) 

C  •••  IF  L  IS  LESS  THAN  M,  U(M)  IS  KNOWN. 

IF(L.EO.M)  U ( M) =GAMMA( M) 

DO  70  I=M, 2,-1 

U ( I- 1 ) aG AMM A( I-1)+R12(I-1)*U(I)/BTA(I-1) 

70  CONTINUE 

RETURN 
END 
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COMPLEX  FUNCTION  HNKL(X) 


•  HANKEL  FUNCTION  H0(1)  -  POLYNOMIAL  APPROXIMATION  • 

•  HANDBOOK  OF  MATH  FUNCTIONS  -  N.B.S.  -  NOV  1967  • 


REAL  JO 

DATA  PI/3. 141592654/ 

IFU.GT.3.)  GO  TO  10 

•••  (-3.0.LE.X.LE.3-0) 

YsX#X/ 9 . 0 

J0=1 .*Y» (-2.2 199997*Y»(*1 .  265  6208*Y'»  ( -0  . 3 163866  *Y»(*0.01 4 1479 
C  ♦ Y* < -0. 00391 «1*Y»<  t.0. 0002  100)  )  )  )  )  ) 


•••  (0.0.LT.X.LE.3.0) 

TO*  2 .  CLOG  (  0 . 5  *X )  »J0/PI*0. 36716691 
C  *Y»( *0.60559366*Y«<-0.79350389*Y»(*0. 25300 1 17*Y«( -0.0426 121 4 
C  ♦ Y»< *0.004279 16* Y«(-0. 00029846) )) ) > ) 

HNKLaCMPLX( JO.YO) 

RETURN 

(3. O.LE.X.LT. INFINITY) 

10  Y*3.0/X 

FO*  0. 797 88 456*Y»(-0. 00000077 *Y» < -0. 00552740+Y»(-0.0n0095 12 
C  *Y»(*0.001 372 3T*Y»( -0.0007 23 05*Y»( *0.00011H7 6) ) ) ) ) ) 

T0*X-0.78539816*Y#(-0. 01166397 ♦Y,(-0.00003951*Y#(*0. 00262571 
C  ♦ Y«(-0. 00051 125*Y»( -0.00029 3 3 3*Y*( *0.0001  3558 ) ) )  ))  ; 

HNKL*FO»CEXP( CMPLX(O.O.TO) )/SQRT(X) 

RETURN 

END 
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SUBROUTINE  S VP< N LYR , ZLYR , RHO , BETA . IXSVP , NSVP , ZSVP , CSVP) 

C 

C  •••  SOUND  VELOCITY  PROFILE  SUBROUTINE 

C  »»•  CALLING  PROGRAM  SUPPLIES:  NOTHING 

C  »»•  SVP  SUBROUTINE  RETURNS: 

C  NLYR  -  NUMBER  OF  LAYERS.  L  YER  1  IS  WATER.  OTHERS  ARE  SEDIMENT 

C  ZLYR  -  ARRAY  -  DEPTH  OF  EACH  LAYER.  FIRST  IS  DEPTH  OF  WATER. 

C  RHO  -  ARRAY  -  DENSITY  OF  EACH  LAYER.  GRAMS/CUBIC  CM 

C  BETA  -  ARRAY  -  ATTENUATION  IN  EACH  LAYER.  DB/WA VELENGTH 

C  IXSVP-  ARRAY  -  CONTAINS  POINTERS.  POINTS  TO  LAST  VALUE  OF  SVP 

C  IN  CORRESPONDING  LAYER.  SVP  IS  STORED  IN  ARRAYS  ZSVP 

C  AND  CSVP.  IXSVP(I)  POINTS  TO  LAST  SVP  POINT  IN  WATER. 

C  IXSVP(NLYR)  POINTS  TO  LAST  SVP  POINT  IN  BOTTOM-MOST  LAYER. 

C  NSVP  -  NUMBER  OF  POINTS  IN  ZSVP  AND  CSVP.  ZSVP  AND  CSVP 

C  CONTAIN  THE  PROFILES  FOR  ALL  LAYERS. 

C  ZSVP  -  ARRAY  -  SVP  DEPTHS  -  METERS 

C  CSVP  -  ARRAY  -  SOUND  SPEED  -  METERS/SEC 

C 

c 

PARAMETER  NIUs1.NPUs6 

DIMENSION  ZLYRC1 ) . RHO< 1 ) , BETA( 1 ), IXSVP (1) , ZSVPC 1 ) .CSVP( 1 ) 

C 

H  S  V  P 1 0 

C  •••  READ  NUMBER  OF  LAYERS 

READ( NIU .» . END= 100)  NLYR 

C  •••  FIRST  LAYER  IS  WATER.  OTHERS  ARE  SEDIMENT. 

DO  55  I  s  1  , NLYR 

C  •••  READ  DEPTH  OF  LAYER.  DENSITY  AND  ATTENUATION. 

R  EAD ( NIU . ENDs 100)  ZLYR( I)  . RHO< I)  . BET  A ( I ) 

C  READ  PROFILE. 

50  READ( NIU ENDs 100)  ZV.CV 

NSVPsNSVP*! 

ZSVP<  NSVP) sZV 
CSVP( NSVP) sCV 

IF( Z V . LT . ZLYR ( I) )  GO  TO  50 
IFC ZV.GT.ZLYRf I) )  GO  TO  100 
IXSVP( I) sNSVP 
55  CONTINUE 

RETURN 
C 

C  •••  ERROR  EXIT 

100  NS VPs  0 

RETURN 
END 
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SUBROUTINE  SFIELD! FRQ. CO , ZS . N, DZ , U) 

C 

C  •••  GAUSSIAN  STARTING  FIELD  -  SEE  NORDA  TECH  NOTE  12  BY  H.K. BROCK 

C  •••  SFIELD  IS  CALLED  IF  INPUT  PARAMETER  ISF  =  0 

C  •■••••••• ••••••••••••••••••••••••••••• 

c 

C  •••  CALLING  ROUTINE  SUPPLIES: 

C  FRQ  -  FREQUENCY  IN  HZ 

C  CO  -  REFERENCE  SOUND  SPEED  -  METERS/SEC 

C  ZS  -  DEPTH  OF  SOURCE  IN  METERS. 

C  N  NUMBER  OF  POINTS  IN  ARRAY  U 

C  DZ  -  DEPTH  INCREMENT  -  METERS 

C  SFIELD  SUBROUTINE  SUPPLIES: 

C  U  COMPLEX  STARTING  FIELD 

c 

COMPLEX  U( 1 ) 

DATA  PI/3. 14 15926535/ 

THE  FIELD  IS  DEFINED  AS  A  GAUSSIAN  BEAM  AT  RANGE  =  0. 

LOCAL  VARIABLES  -  GA  GAUSSIAN  AMPLITUDE 

XK0=2.0»PI»FRQ/C0 

GW=2 . O/XKO 

GAsSQRT( GW) /GW 

DO  10  1=1. N 

2M=I»DZ 

PR=GAUSS(GA,ZM,ZS.GW)-GAUSS(GA,-ZM,ZS,GW) 

U(I)=CMPLX(PR.O.O) 

0  CONTINUE 
RETURN 
END 

FUNCTION  GAUSS(GA.Z.GD.GW) 

INPUT  -  GA  GAUSSIAN  AMPLITUDE 

OUTPUT  -  GAUSS  =  GA  •  EXP(-((Z  -  GD)  /  GW)»»2) 

TEMPORARY  VARIABLE  -  V 
V=( Z-GD)/GW 
V=-( V«V) 

GAUSS=GA»EXP( V) 

RETURN 

END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 

100 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE  USVP 


•••  USER  SOUND  VELOCITY  PROFILE  SUBROUTINE 

SUBROUTINE  USVP  IS  CALLED  EACH  DR  IN  RANGE  AS  LONG  AS 

KSVP  IS  NOT  ZERO.  KSVP  MAY  BE  USED  BY  USER  TO  TRANSFER  CONTROL 

IN  THIS  SUBROUTINE.  USER  INSERTS  LOGIC  TO  CLEAR  KSVP 

WHEN  USVP  IS  NO  LONGER  NEEDED.  IF  KSVP  NOT  CLEARED  BY  USER, 

USVP  IS  CALLED  EACH  STEP  IN  RANGE  UNTIL  RA  s  NEXT  RSVP. 


•••  USVP  SUBROUTINE  RETURNS: 

NLYR  -  NUMBER  OF  LAYERS.  LAYER  1  IS  WATER.  OTHERS  ARE  SEDIMENT 
ZLYR  -  ARRAY  -  DEPTH  OF  EACH  LAYER.  FIRST  IS  DEPTH  OF  WATER. 
RHO  -  ARRAY  -  DENSITY  OF  EACH  LAYER.  GRAMS/CUBIC  CM 
BETA  -  ARRAY  -  ATTENUATION  IN  EACH  LAYER.  DB/W A VELENGTH 
IXSVP  -  ARRAY  -  CONTAINS  POINTERS.  POINTS  TO  LAST  VALUE  OF  SVP 
IN  CORRESPONDING  LAYER.  SVP  IS  STORED  IN  ARRAYS  ZSVP 
AND  CSVP.  IXSVP(I)  POINTS  TO  LAST  SVP  POINT  IN  WATER. 

NS VP  -  NUMBER  OF  POINTS  IN  ZSVP  AND  CSVP.  ZSVP  AND  CSVP 
CONTAIN  THE  PROFILES  FOR  ALL  LAYERS. 

ZSVP  -  ARRAY  -  SVP  DEPTHS  -  METERS 
CSVP  -  ARRAY  -  SOUND  SPEED  -  METERS/SEC 
KSVP  -  AS  DESCRIBED  ABOVE. 


PARAMETER  MXLYR* 101 . MXN*  1 0000 . MXSVP* 1 0 1 , MXTRK* 1 0 1 . NIU=  1  . 

C  NOU*  2 , N  PU  =  6 

COMPLEX  ACOFX, ACOFY, BCOF.BOTX, 80TY . BTA . HNK , HNKL . SURX . SUR Y . TEMP , 

C  U.X.Y 

COMMON  /IFDCOM/ACOFX, ACOFY, ALPHA. BCOF , BETA( MXLYR ) .BOTX.BOTY. 

BTA( MXN )  , CO. CSVP( MXSVP) , DR , DR  1 . DZ , FR Q . IHNK . ISF , ITYPEB , 
ITYPES, IXSVP( MXLYR) , KS VP , N , N 1 , NLYR . NSVP . NWSVP , R 1 2 ( MXN ) , R A , 
RHO(MXLYR)  , RSVP , SUR X . SUR Y , THETA , TRACKC MXTR K  .  2 )  ,U(MXN)  , 

X ( MXN ) , XKO , Y ( MXN ) , Z A , Z LYR ( MXL YR ) , Z P , ZS . ZS VP ( MXS VP) 

DATA  P 1/3.  I'll  5  9265  H/  ,  DEG/5  7. 2957  8/ 

GO  TO  (  100.200,300.1100)  .KSVP 

NSVPsO 

RETURN 

CONTINUE 

IF  KSVP* 1 ,  CONTROL  IS  TRANSFERRED  HERE.  USER  LOADS 
NLYR.ZLYR(I) ,RHO(I)  , BET A{ I )  ,  AND  IXSVP(I)  WHERE  1*1, NLYR. 

USER  ALSO  LOADS  NSVP , ZS VP( I ) .  AND  CSVP(I)  WHERE  1=1, NSVP. 

KSVP  MAY  BE  ALTERED  DEPENDING  ON  USER  LOGIC. 

•••  USER  SUPPLIES  SVP 

»*•  SETUP  FOR  ONE  LAYER  WITH  FOUR  SVP  POINTS  SHOWN  BELOW 
NLYR* 1 

ZLYR ( 1)s . 

RHO( 1 ) *  1 . 0 

BET  A ( 1 ) * . 

NSVPsfl 


I 


A- 31 


o-ero  oiAiO  n  fu  n  ooooooo 


TR  6659 


C  ZSVP(1)s . 

CSVP( 1 )s . 

ZSVP( 2) s . 

CSVP( 2 ) = . 

ZSVP( 3) . 

CSVP( 3 ) - . 

ZSVPU)  = . 

CSVP( «) = . 

IXSVPC 1 ) sNSVP 
RETURN 

00  CONTINUE 

USER  INSERTS  CODE  HERE  IF  DESIRED 
RETURN 

00  CONTINUE 

•••  USER  INSERTS  CODE  HERE  IF  DESIRED 
RETURN 

00  CONTINUE 

USER  INSERTS  CODE  HERE  IF  DESIRED 
RETURN 
END 
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SUBROUTINE  UFIELD 

c  •*•••••••••••••••••••••••••••••••••••*••*•••••*•••••••••••••••• 

C  USER  STARTING  FIELD 

C  •••  USER  WRITES  THIS  SUBROUTINE  IF  GAUSSIAN  FIELD  NOT  DESIRED 

C  •••  UFIELD  IS  CALLED  IF  INPUT  PARAMETER  ISF  IS  NOT  ZERO 

C  •••  UFIELD  SUBROUTINE  SUPPLIES: 

C  U  COMPLEX  STARTING  FIELD 

C  •••••••••••••••••••••••••••••••••••#•••••••••••••••••••»•••#••• 

PARAMETER  MXLYRs  1  0  1  .  MXNs  ’TOOO  ,  MXS VP  =  1  0  1  , MXTRKs  101  . NIUs 1 , 

C  NOU*  2 , NPUs  6 

COMPLEX  AIOFX.ACOFY.BCOF.&OTX.BOTY.BTA.HNK.HNKL, SURX , SU R Y . TEM P , 

C  U.X. Y 

COMMON  /IFDCOM/ACOFX.ACOFY  ALPHA , BCOF , 8ETA( MXLYR ) , BOTX, BOTY, 

BTA( MXN) . CO, CSVP( MXS VP) . DR . DR  1 . DZ . FRO, IHNK, ISF, ITYPEB . 

I  TYPES, IXS VP ( MXLYR)  , KSVP. N. N1 . NL YR . NS VP , NWS VP . R 1 2( MXN)  , RA , 
RHO( MXLYR) , RS V P , SUR X . S UR Y . THET A , TR ACK( MXTR K , 2 ) ,U(MXN)  , 
X(MXN) ,XKO,Y( MXN) ,ZA,ZLYR( MXLYR ) ,ZP,ZS.ZSVP( MXS VP) 

DATA  PI/3. 14159265N/. DEG 57 8/ 

C 

C  STARTING  FIELD  GENERATED  BY  USER 

DO  10  Isl.N 
Z  I  s  I  *  D  Z 

C  U(I)= . 

10  CONTINUE 

RETURN 
END 
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SUBROUTINE  BOON 

C  » 

C  USER  PREPARED  BOTTOM  CONDITION  SUBROUTINE 

C  BOON  IS  CALLED  IF  INPUT  PARAMETER  ITYPEB  =  1 

C  SEE  MAIN  PROGRAM  FOR  DEFINITIONS 

C 

C  •••  SUBROUTINE  RETURNS: 

C  BOTT.BOTX 

c  •••■■•••••••»■••••••■••••••••••••••••••••••••••••••••••••«••••••• 

c 

PARAMETER  MXLYRs 1 0 1 , MXNs 1 OOOO . MXS VPs  1 0 1 , MXTRKslOI , NIUsI , 

C  NOUs  2 , N  PUs  6 

COMPLEX  ACOFX, ACOFY, BCOF, BOTX. BOTY , BTA , HNK , HNKL , SU R X , SU R Y . TEMP . 

C  U.X.Y 

COMMON  /IFDCOM/ ACOFX, ACOFY, ALPHA, BCOF , BETA( MXLYR ) , BOTX , BOTY , 
BTA(MXN) , CO , CSVP( MXSVP) .  DR  ,  DR  1  ,  DZ  .  FR  Q ,  I HN  K ,  ISF  ,  ITY PEB  , 

I TYPES. IXSVP< MXLYR) , KS VP . N , N 1 , NL YR , NS VP , NWS VP , R 1 2 ( MXN ) , RA . 
RHO(  MXLYR) , RS V P , SUR X . SUR Y , THETA , T R ACK( MXTR K . 2 ) , U(MXN)  , 

X ( MXN ) , XKO.Y(MXN) , ZA , ZLYR( MXLYR) , ZP . ZS , ZS VP( MXS VP ) 

DATA  PI/ 3.  1<*  15 9265**/  .  DEG/ 5^  -29578/ 

C 

IF ( THETA )  50,100,150 
C 

C  •••  THETA  LESS  THAN  0.0,  BOTTOM  SLOPES  UP. 

50  CONTINUE 

BOTYsU(N) 

C  BOTX  s . 

RETURN 

C 

C  •••  THETA  s  0.0.  BOTTOM  IS  FLAT. 

100  CONTINUE 

BOTYs  U ( N ) 

C  aOTXs . 

RETURN 

C 

C  •••  THETA  GREATER  THAN  0.0,  BOTTOM  SLOPES  DOWN. 

150  CONTINUE 

C  BCTYs . 

C  BOTX  s . 

RETURN 

END 
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SUBROUTINE  SCON 

llttHlltllllHIHHHIHIMIIItMlltlllHHMMIIItmi 

•••  SURFACE  CONDITION  SUBROUTINE 

IF  ITYPE3  ■  0.  SCON  SETS  SURY  AND  SURX  >  0.0. 

IN  ITT DCS  NOT  0.  THE  USER  MUST  SUPPLY  SURY  AND  SURX. 

Sit  MAIN  PROGRAM  POR  DEFINITIONS 


PARAMETER  NXLTIa 101, MXN.1 0000, MXSVPs 1 0 1. MXTRKs 1 0 1, NIUs 1 , 

C  NOUa2«  NPU»6 

COMPLEX  ACOFX . ACOFY, BCOF , BOTX , BOTY, BTA , HRK, HNKL , SURX . SURY . TEMP , 

C  O.X.Y 

COMMON  /IFDCOM/ ACOFX. ACOFY, ALPHA, BCOF, BETA(MXLYR) .BOTX, BOTY. 

C  STA(MXN) .CO.CSVP(MXSVP) , DR, DB 1 , D2.FR0, IHNK. ISF . ITYPEB , 

C  ITTPKS, IXSVP( MXLYR) , XSVP, N, N 1 , NLYR , NS VP, NWSVP, R12( MXN) ,RA, 

C  RNO< MXLYR). RSVP, SURX. SURY, THETA. TRACKC MXTRK, 2)  ,U(MXN)  . 

C  X(MXN) , XXO, Y( MXN) , ZA ,ZLYR( MXLYR) , ZP, ZS,ZSVP( MXSVP) 

DATA  PI/3. 1*159265*/. DEG/57. 29578/ 

C 

IF( I TYPES. RE. 0)  CO  TO  100 
C 

C  •••  PRESSURE  RELEASE  SURFACE 

SOSTaO.O 
SURXaO.O 
RETURN 

•••  USER  SURFACE  CONDITION 
100  CONTINUE 

C  SURY . 

C  SURX . 

RETURN 

INO 
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Appendix  B 

PLOT  PROGRAM  COMPUTER  LISTING 


B-l/B-2 
Reverse  Blank 
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***  PLOT  PROGRAM  POR  I PD  MODEL. 

***  PROGRAM  PLOTS  PROPAGATION  LOSS  VS  RANGE  ON  CALCOMP  PLOTTER  - 
***  MODEL  1039. 

**<*  SEE  MAIN  PROGRAM  IPD  FOR  DEFINITIONS  OP  IFD  VARIABLES. 


**•  INPOT 


INPUT  UNIT  NUMBER  -  NIU 
INPUT  PILE  NAME  ■  PLTIFD.IN 
CONTENTS:  CARD  IMAGES  IN  FREE  FORMAT 


CARO  1 
CARD  N 


Z1,Z2,Z3,Y1,Y2,Y3,YL,X1,X2,X3,XL,PACT,XAVG 


♦**  QUICK  REFERENCE  AND  NOTES  FOR  CARD  INPUT 


21  «  FIRST  RECEIVER  DEPTH  TO  PLOT  -  METERS 
IP  (Z1.LE.0.0)  PROGRAM  TERMINATES 
Z2  ■  LAST  RECEIVER  DEPTH  TO  PLOT  -  METERS 
Z3  «  RECEIVER  DEPTH  INCREMENT  -  METERS 
Y1  •  LABEL  OF  Y-AXIS  AT  ORIGIN  IN  DB 
Y2  -  LABEL  AT  TOP  OF  Y-AXIS  IN  DB 
Y3  -  INCREMENT  OF  Y-AXIS  LABELS  TN  DR 
YL  -  LENGTH  OF  Y-AXIS  IN  INCHES 
XT  -  LAPEL  OP  X-AXTS  AT  ORTOTN  IN  «fTL9MSTf\RS 
X2  -  LABEL  AT  RIGHT  OF  X-AXIS  IN  *TLC*FTERS 
X3  -  INCREMENT  OF  X-AXIS  LABELS  IN  KILOMETERS 
XL  -  LENGTH  OF  X-AXIS  IN  INCHES 

FACT  -  SCALE  OF  PLOT:  1.0  »  FULL  SIZE  ;  .5  -  1/2  SIZE  »  BTC. 
XAVG  •  RANGE  OVER  WHICH  TO  COMPUTE  RUNNING  AVERAGE  IN  METERS 
IF  XAVG  ■  0,  ALL  POINTS  ARE  PLOTTED 


PARAMETER  MAXP-5000 

PARAMETER  MXLYR- 1 0 1 , MXN-1 0000 , NIU-1 , NOU-2 , NPU-S , PLTU-3 
COMPLEX  HNK,HNXL,CTBMP,U(MXN) 

DIMENSION  BETA (MXLYR) ,RHO (MXLYR) , ZLYR (MXLYR) ,IBUF(2000) 
DIMENSION  P(MAXP) ,R(MAXP) 

DATA  PI/3. 141592*554/ ,  DEG/57. 29578/ 

DATA  CNVKM/1000.0/ 

IPRNT-0 

C  ***  ASSIGN  IFD  OUTPUT  FILE 

CALL  ASSIGN <NOU,' IFD. OUT' ) 

C  ***  ASSIGN  PLOT  PARAMETER  INPUT  FILE 
CALL  ASSIGN (NIU, 'PLTIFD. IN' ) 

CALL  PLOTS<IBUF,2000,PLTU) 

CALL  PLOT  (0.0,0. 5,-3) 

C  ***  READ  PLOT  PARAMETERS 

100  READ(NIU,*,ENO-510)  Z1,Z2,Z3,Y1,Y2,Y3,YL,X1 ,X2 ,X3 , XL, FACT, XAVG 
IP(Zl.LE.O.O)  GO  TO  510 
I F ( FACT • LE . 0 . 0 )  FACT-1.0 
CALL  FACTOR (FACT) 
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YINC-(Y2-Y1)/YL 
IF(Z3.E0.0.0)  Z3-1.0 

C  ***  GENERATE  PLOT  FOR  RECEIVER  OEPTH 
DO  350  ZR-Z1,Z2,Z3 
ZRR-ZR 
I  PEN- 3 
IX-1 

DX-(X2-X1)/XL 

CALL  AXIS2(0. ,0. ,* RANGE  (KM) ( ,-10,XL,C.,Xl,X3,X2) 

CALL  AXIS2 (XL,0. , '  • ,-l ,YL,90. ,Y1 ,Y3 ,Y2) 

CALL  AXIS2 (0 . , 0.  ,  ' PROPLOSS  (DB) ' ,+13 . YL,90. ,Y1 ,Y3 ,Y2) 

CALL  PLOT(0.0,YL,3) 

CALL  PLOT (XL, YL, 2) 

C  ***  READ  INITIAL  IPD  PARAMETERS 

REWIND (NOU) 

READ(NOU)  FRQ#ZS,C0 , ISP,R0 ,ZP,N, IRNK, ITYPBB , ITYPES,RWAX,DR»WDR,D7 , 
CNLYR , ZLYR , RHO , BETA 
IF(XAVG.LT.WDR)  XAVG-WDR 
L-0 

115  CONTINUE 
RAVG-0.0 
PLAVG-0 . 0 

C  ***  READ  SOLUTION  FIELD 

120  READ (NOU , END— 170) NN,RA ,WDZ , (U (I ) , I— 1 ,NN) 

IP(RA.LE.O.O)  GO  TO  120 
HNK-HNKL (2 . 0*PI*FRO*RA/C0) 

1-0 

INTERP— 0 
I -ZRR/WDZ 

IF(I.GT.NN)  GO  TO  115 
IF(I.GE.l)  GO  TO  130 
1-1 

ZRR-WDZ 

130  IP(I*WD2.NE.ZRR)  INTERP-1 
Y-CABS (U (I)  ) 

IF(IHNK.NE.O)  Y-CABS (U (I )*HNK) 

I F ( INTERP. EQ. 1 )  CTEMP-U (I ) ♦ (U ( I+l ) -U ( I ) ) * ( ZRR-I*WDZ) /WDZ 
IF ( INTERP. EQ. 1. AND* IRNR, ME. 0)  Y-CABS (CTEMP*HNK) 

IF(Y.LE.O.O)  GO  TO  120 
L-L+l 
P  (L)  -Y 
R(L)-RA 
GO  TO  120 
170  CONTINUE 
K-0 

200  K-K+l 
RAVG-0 
PLAVG-0 
NAVG-0 

DO  210  «J-K ,  L 

IF(R(J)-R(K) .CE.XAVG)  GO  TO  220 
RAVG-RAVG+R(J) 

PLAVO-PLAVG+P (J) 

MAVG-WAVG*1 
210  CONTINUE 
220  CONTINUE 

IF(MAVG.EQ.O)  GO  TO  250 
BIAS-0.0 
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RA-RAVG/NAVG 

IF(IHNK.EQ.O.ANO.RA.GT.O.O)  BIAS-lO.O*ALOG10 (RA) 

Y-PLAVG/NAVG 
IF(Y.LE.O.O)  GO  TO  200 
Y--2Q.0-ALCG10  (Y)+BIAS 
TEMP-RA/ l OOO.fi 

IF(IPRNT.EQ.l)  WRITE (MPU,*)  TEMP,Y 

IF (RA/CNVKM.GT.X2-XAVG/CNVKM)  GO  TO  250 

IF(RA/CNVKM.LT.X1)  GO  TO  200 

X* (RA/CNVKM-X1 ) /OX 

Y-(Y-Y1)/YINC 

IF(Y.LT.O.O)  Y«0.0 

TF(Y.GT.YL)  Y-YL 

CALL  PLOT(X,Y,IPEN) 

I  PEN-2 
GO  TO  200 
250  CONTINUE 

CALL  BLOCK(XL,YL,FRO,ZS,CO,ISF,RO,ZO,V,IHNK,ITYPEB,ITYPES, 

CRMAX , DR , WOR, WDZ , DZ , ZRR ,NLYR ,  ZLYR , BETA , RHO , XAVG ) 

CALL  PL0T(XL+2. 0,0. 0,-3) 

350  CONTINUE 
GO  TO  100 
510  CONTINUE 

CALL  PLOT (0.0,-0.5,999) 

STOP 

END 

SUBROUTINE  BLOCK (XL , YL , FRQ, ZS ,C0 , ISP , RO , ZO , N, IHNK , ITYPEB , ITYPES , 
CRMAX , DR ,WDR , WOZ , OZ , ZRR , NLYR , ZLYR , BETA , RHO , XAVG ) 

DIMENSION  Z LYR ( 1 ) , BETA ( 1 ) , RHO ( 1 ) 

MC-30  !  MAX  CHAR  IN  STRING 

HT-.l 

DY-1.5*HT 

XBLK«5.0*HT 

YBLK-YL+ (21+NLYR) *DY 

CALL  SYMBOL(XBLX,YBLX,HT,*TFD  SOLUTION* ,0. ,12) 

YBLK-YBLK-DY 

CALL  SYMBOL (XBLK,YBLK,HT, *  INITIAL  PARAMETERS* ,0.0,18) 

YB  LK  — ' YB  LK — DY 

CALL  SYMBOL (XBLK,YBLK,HT,* PRO  ■  ',0.,S) 

CALL  NUMBER (999. ,YBLK,HT,FRO,0. ,1) 

CALL  SYMBOL (999. ,YBLK,HT,'  HZ*,0.,3) 

YBLK-YBLK-DY 

CALL  SYM80L(XRLK,YBLK,HT,'ZS  ■  *,0.,5) 

CALL  NUMBBR(999.,YBLK,HT,ZS,0.,1) 

CALL  SYMBOL (999. ,YBLX,HT,*  M*,0.,2) 

YBLK-YBLK-DY 

CALL  SYMBOL (XBLK , YB  LK , HT , ' CO  -  *,0.,5) 

CALL  NUMBER(999. ,Y8LK,HT,C0 ,0. ,1) 

CALL  SYMBOL (999. ,YBLK,HT,*  M/SEC*, 0., 5) 

YBLK-YBLK-DY 

CALL  SYMBOL(X8LK,YBLK,HT,*RO  -  ’,0.,5) 

CALL  NUM8BR(999.,YSLK,HT,R0,0.,1) 

CALL  SYMBOL (999. ,YBLK,HT,*  M*,0.,2) 

YBLK-YBLK-DY 

CALL  SYMBOL (XBLK, YBLK,HT, ' ZO  •  ',0.,5) 

CALL  NUMBER (999. ,YRLX,HT,Zfl,0. ,1) 

CALL  SYMBOL (999. ,YBLK,HT,'  M*,0.,2) 
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Y8LX-YBLK-0Y 

CALL  SYMBOL (XB  LX  ,  YB  LK ,  HT , '  N  -  • ,0. ,5) 

FN-N 

CALL  NUM8ER(999. ,YRLK,HT,FN,0. ,-\) 

YBLK-YBLX-DY 

CALL  SYMBOL  (XBLX, YBLK,  HT, '  OR  -  5) 

CALL  NUMBER (99°. , YBLK,HT,DR,0. ,1) 

CALL  SYMBOL (999. , YBLK, HT, '  M',0.,2) 

YQLK*YBLK-OY 

CALL  SYM80L(XBLK,YBLK,  HT, 1  WDR  •  ',(>., 9) 

CALL  NUMBER (999. , YBLK, HT, WDR, fl. ,1) 

CALL  SYMBOL (999. , YBLK,  Iff,'  M',0.,2) 

YBLK-YBLK-DY 

CALL  SYMBOL (XBLK,YBLK,HT,'RMAX  ■  *,0.,7) 

CALL  NUMBER (999. , YBLK, HT,RMAX,0. ,1) 

CALL  SYMBOL (999. ,YBLX,HT,'  M',0.,2) 

YQ  LK» YB LX *DY 

CALL  SYMB0L(XBLK,Y8LK,HT,’DZ  -  *,0.,5) 

CALL  NUMBER(999.  ,YBLK,HT,DZ,0.  ,2) 

CALL  SYMBOL (999. ,YB LX, HT,*  M',0.,2) 

YBLK-YBLK-DY 

CALL  SYMBOL  (XBLK.YB  LX,  HT,'WDZ  *  *,0.,<) 

CALL  NUMBER(999. ,YBLK,HT,WDZ ,0. ,2) 

CALL  SYMBOL (999. ,YBLK,HT,'  M’,0.,2) 

YBLK-YBLX-DY 

CALL  SYMBOL (XBLK, YBLK, HT,' ISP  -  ’,0.,*) 

FPN-ISF 

CALL  NUMBER (999. ,YBLX, HT,FPN,0. ,-l) 

YBLK-YBLX-DY 

CALL  SYMBOL (XBLK.YB LX, HT,' I HNK  -  ' ,0. ,7) 

FPN-IHNK 

CALL  NUMBER (999. ,YBLK,HT,FPN,0. ,-I) 

YBLX-YBLK-DY 

CALL  SYMBOL (XBLK.YB LX, HT.'ITYPES  -  *,0.,9I 
PPN-ITYPES 

CALL  NUMBER (999 . , YBLK , HT , FPN ,0. ,-l) 

YBLK-YBLX-DY 

CALL  SYMBOL(XBLX, YBLK, HT, * TTYPEB  -  ’,0.,9) 

PPN-ITYPEB 

CALL  NUMBER (999., YBLK, HT,PPN,0.,-1) 

YBLK-YBLX-DY 

CALL  SYMBOL (XBLK, YBLK, HT.'LYR  DEPTH (M)  RHO  BETA(DB/WL) * ,0. ,12) 

DO  son  I-l.NLYR 

YBLK-YBLK-DY 

FPN-I 

CALL  NUMBER (XBLK, YBLK, HT, FPN, n. , -I) 

CALL  NUMBER (XBLK+5. 0*HT, YBLK, HT,ZLYR(I> ,0.,1) 

CALL  NUMBER (XBLX+1 4. 0*HT, YBLK, HT, RHO (T) ,0.,2) 

IF(BETA(I) .GE. 9.0) CALL  NUMBER (XBLK+71 .0*HT, YBLK, HT.BETA (I) ,0.,3) 
IF(BETA(T) . LT.O,)CALL  SYMBOL (XBLK+21 . n*HT, YBLK, HT, 'COMPUTED' ,0. ,P) 
900  CONTINUE 

YBLK-YBLK-DY 

CALL  SYHBOL(XBLK, YBLK, HT, ' AVG  -  ’,(*.,7) 

CALL  NUMBER { 999 ., YBLK , HT , XAVG , 0 . , - 1 ) 

CALL  SYMBOL (999. , YBLK, HT, '  M',0.,2) 

YBLK-YBLK-DY 

CALL  SYMB0L(X8LX, YBLK, HT, 'RECEIVER  DEPTH  -  ',0.,17) 

CALL  NUMBER ( «9  9 . , YBLK , HT , ZRR , 0 . , 1 ) 

CALL  SYMB0L(999.  .YBLK.HT, '”  M»,0.,2) 

RETURN 

END 
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